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Preface 


The  purpose  of  this  study  was  to  validate  the  three-body 
perturbation  theory  described  by  Dr.  William  Wiesel,  in  his  paper 
entitled,  "Perturbation  Theory  in  the  Vicinity  of  a  Periodic 
Orbit  by  Repeated  Linear  Transformations".  Mastering  these 
techniques  would  be  a  necessary  first  step  to  navigating  a  highly 
perturbed  region  of  space,  such  as  that  surrounding  the  Martian 
moon  Phobos.  Any  further  developments  in  the  understanding  of 
time-periodic  systems  in  general,  could  reap  great  rewards  in 
several  non-astronautical  areas  as  well.  The  dynamics  of 
helicopter  blades  is  a  good  example. 

This  entire  study  would  not  have  been  possible  if  not  for 
the  extensive  help  I  received  from  Dr.  Wiesel  himself.  The 
process  began  with  many  hours  of  tutoring  on  the  major  technical 
points  of  the  thesis.  Each  area  discussed  eventually  became  one 
of  the  six  different  programs  used  during  the  course  of  the 
study.  I  was  given  free  reign  of  several  different  programs, 
subroutines,  and  hard  to  find  text  books.  If  Dr.  Wiesel  didn't 
already  have  a  similar  program  to  be  modified,  he  would  help  me 
lay  out  an  algorithm  that  I  could  later  encode.  Dr.  Wiesel  also 
devised  several  analytical  tests  that  could  be  made  to  ensure 
program  accuracy. 

For  all  of  the  reasons  stated  above,  as  well  as  the  thrill 
of  showing  someone  that  their  ideas  actually  work,  I  consider 
this  thesis  an  overwhelming  success.  Yet  there  are  still  a  great 
deal  of  unanswered  questions,  and  program  modifications  to  make. 
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I  would  also  like  to  take  this  opportunity  to  thank  my  wife 
Beth,  for  her  support  and  understanding  during  this  last  year  and 
a  half.  Without  her,  the  answer  to  the  question,  why  bother  ?, 
would  not  have  been  so  apparent. 

Finally,  I  hereby  dedicate  this  entire  work  to  my  very 
special  other  Mom,  Lynne,  whose  passing  has  opened  my  eyes  to 
what  is  truly  important  in  this  world. 

David  A.  Ross 
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ABSTRACT 

A  perturbation  theory  for  restricted  three-body  orbits, 
using  a  periodic  trajectory  as  a  reference  solution,  is 
investigated.  The  nearly-periodic  equations  of  motion  are 
derived  by  analogy  to  a  linearization  about  an  equilibrium  point 
In  this  case,  the  linearization  produces  a  set  of  time-periodic 
equations  of  motion  that,  according  to  Floquet,  are  completely 
solved  by  a  periodic  trajectory. 

The  four-dimensional  phase  space  of  the  restricted  three- 
body  problem  is  first  surveyed  for  regions  of  periodic  motion, 
via  the  surface  of  section  phase  plot.  Upon  extraction  of  a 
periodic  orbit,  nearly-periodic  orbits  are  integrated.  The 
integrated  state  vector  is  routinely  sampled,  and  then  twice 
transformed  into  a  new  set  of  variables.  The  first  translates 
the  frame  center  to  the  periodic  trajectory.  The  second,  or 
modal  transformation,  projects  the  coordinates  along  their 
eigenvectors.  The  transformations  are  highly  useful,  since  two 
of  four  new  variables  are  constant  within  a  finite  region- 
surrounding  the  periodic  reference.  Plots  of  the  two  variables 
are  offered  as  an  exact  representation  of  a  nearly-periodic 
trajectory,  while  plots  of  the  constants  over  time,  trace  the 
boundaries  of  the  nearly-periodic  region.  ^ 

After  the  original  Hamiltonian  is  canonically  transformed 
into  the  new  variables,  it  is  expanded  in  a  Taylor's  series. 
Several  of  the  terms  are  either  simplified  or  annihilated 
completely.  The  expansion  is  then  truncated  after  four  terms. 
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leaving  a  readily  differentiable  expression  from  which  to  derive 
the  nearly-periodic  equations  of  motion.  The  expanded 
trajectories  are  then  compared  to  the  exact  ones,  over  a  wide 
range  of  values. 

As  was  expected,  a  significant  region  exists  where  the 
expanded  equations  of  motion  accurately  reproduce  the  concentric 
circular  paths  shown  to  exist  by  the  transformed  case.  As  the 
initial  displacements  from  the  periodic  trajectory  are  increased, 
the  expanded  trajectories  fail  to  accurately  model  the  transients 
observed  in  the  exact  case.  These  exact  case,  orbital 
irregularities  occur  because  the  displacements  from  the  periodic 
orbit  are  no  longer  small  enough  to  be  represented  by  a 
linearization  of  the  dynamics.  The  expanded  trajectories  fail  to 
recreate  this  non-linear  behavior  because  most  of  the  Hamiltonian 
terms  responsible  have  been  truncated.  Therefore,  before  the 
complete  perturbation  solution  can  be  constructed,  the  expanded, 
nearly-periodic  equations  of  motion  should  be  derived  again  using 
more  than  four  terms  of  the  expanded  Hamiltonian. 
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PERTURBATION  THEORY  FOR  RESTRICTED  THREE-BODY  ORBITS 


I.  Introduction 

Before  the  astrodynamics  of  man-made  objects  in  space  can  be 
fully  understood,  one  must  first  comprehend  basic  planetary 
motion.  Thanks  to  Sir  Isaac  Newton  and  his  three  laws  of  motion, 
and  to  Johann  Kepler  for  his  three  laws  of  orbital  motion,  it  can 
be  shown  that  nearly  all  astrodynamical  systems  are  dominated  by 
a  single  conservative  force  known  as  gravity.  In  fact,  the  most 
general  description  of  the  motion  of  a  collection  of  objects  in 
space  is  defined  by  the  n-body  problem. 

In  an  n-body  system,  the  n'**  body  is  acted  upon  by  the  other 
n-1  gravitational  masses  present.  In  this  way,  the  motion  of  any 
mass  in  a  system  affects  and  is  affected  by  every  other  mass  in 
the  system.  The  overwhelming  task  of  representing  each  body  is 
well  illustrated  by  Wiesel. 

Our  own  solar  system  consists  of  one  star,  nine  planets, 
over  fifty  moons,  tens  of  thousands  of  asteroids,  and 
millions  of  comets.  The  description  of  the  motion  of  this 
system  is  clearly  important,  but  an  exact  solution  to  this 
problem  has  not  been  found  in  over  three  hundred  years  of 
study.  (8:33) 

Therefore,  the  use  of  the  exact  n-body  description  of  a  dynamical 
system  is  not  simply  a  nuisance,  it  is  virtually  impossible  to 
implement . 

The  simplest  and  most  drastic  approximation  to  the  n-body 
problem  is  known  appropriately  as  the  two-body  problem.  Here, 
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only  two  point  masses  are  considered  to  exist  within  the 
dynamical  system.  The  primary  masses  are  then  constrained  their 
by  mutual  gravitational  attraction.  Specifically,  both  bodies 
remain  in  a  circular  orbit  about  the  system  center  of  mass  point. 
This  particular  arrangement  is  very  special  in  the  field  of 
astrodynamics  since,  "it  is  the  only  gravitational  problem  for 
which  a  closed-form  solution  has  been  found"  (8:45). 

Additionally,  one  might  suspect  the  accuracy  of  such  a  severely 
truncated  formulation.  Surprisingly  enough,  "most  systems 
encountered  in  orbital  mechanics  are  nearly  perfect  two-body 
problems,  with  only  small  perturbations  from  two-body  motion" 
(9:75)  . 

There  are  situations,  unfortunately,  where  the  small 
perturbation  assumption  of  two-body  perturbation  theory  is 
violated.  Such  a  case  occurs  in  the  vicinity  of  the  Martian  moon 
Phobos.  The  orbit  about  Phobos  is  so  dramatically  effected  by 
the  gravitational  pull  from  Mars,  that  a  simple  two-body 
approximation  can't  simulate  the  true  dynamics  of  the  region. 
According  to  Szebehely, 

Entry  into  celestial  mechanics  and  space  dynamics  can  be 
gained  by  the  study  of  the  problem  of  two  bodies.  To 
penetrate  the  fundamental  problems,  the  number  of 
participating  bodies  must  be  increased  from  two  to  three. 
(5:v) 

Clearly,  in  highly  perturbed  dynamical  systems,  higher  orders  of 
approximation  must  replace  the  tractable  two-body  scenario. 
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II.  Historical  Development 


This  study  is  the  first  of  it's  kind.  An  extensive 
literature  search  of  four  data  bases,  the  AFIT  library,  and  a 
tedious  bout  with  the  Science  Citation  Index,  produced  nothing. 
The  only  sources  of  information  available  on  the  development  of  a 
restricted  three-body  perturbation  theory  using  a  periodic  orbit 
as  the  reference  solution,  were  Dr.  Wiesel's  paper  on  the  subject 
and  Dr.  Wiesel  himself.  As  described  in  his  previous  work.  Dr. 
Wiesel  contends  that. 

It  is  common  for  researchers  working  with  periodic  orbits, 
to  also  solve  the  associated  Floquet  problem  in  order  to 
derive  stability  information  on  the  orbit.  It  is  far  less 
common  to  make  use  of  eigenvectors  of  the  linearized  system, 
or  to  use  a  periodic  orbit  as  a  reference  solution  for 
perturbation  theory.  (7:231) 

While  it  appears  that  this  study  is  similar  to  others  in  it's 
content,  it  is  quite  original  in  it's  purpose.  It  is  unique  to 
use  the  eigenvectors  and  Poincar6  exponents  of  the  periodic 
trajectory,  to  canonically  transform  the  generic  equations  of 
motion  into  nearly-periodic  ones.  By  analogy  to  classical  two- 
body  perturbation  theory  the  periodic  trajectory  serves  as  the 
reference  solution,  while  displacements  from  this  reference  are 
treated  as  small  perturbations. 
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III.  Theory 


The  Restricted  Three-Body  Problem 

The  complete  three-body  formulation  at  first  glance,  might 
not  appear  to  be  much  more  difficult  to  solye  than  the  two-body 
problem.  While  the  inclusion  of  another  grayitational  mass  into 
a  two-body  system  would  further  complicate  the  dynamics,  one 
would  still  hope  to  find  at  least  a  partial  solution.  In 
reality,  however,  the  problem  of  three-bodies  is  completely 
unsolvable  without  imposing  the  restrictions  first  defined  by 
Leonard  Euler  in  1772. 

In  the  restricted  three-body  problem,  it  is  assumed  that  two 
of  the  bodies  are  tremendously  more  massive  than  the  third.  The 
motion  of  the  third  body  is  then  governed  by  the  gravitational 
pull  of  the  two  primary  bodies.  Conversely,  the  motion  of  the 
primaries  is  unaffected  by  the  third  body,  and  is  completely 
described  by  the  two-body  solution.  In  this  way,  the  restricted 
problem  can  be  considered  a  one-body  problem,  because  only  the 
equations  of  motion  of  the  third  body  are  of  interest. 

Before  deriving  the  equations  of  motion,  non-dimensional 
variable  definitions  for  mass,  length,  and  time  must  be 
introduced.  First,  the  non-dimensional  masses  mi,  mj,  and  m3  are 
defined  by 


% 


Ml  M2 


AT, 


Mi^M2 


•0 


(1) 


V7here  Mi  and  M2  are  the  masses  of  the  primary  bodies,  and  M^  is 
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the  mass  of  the  third  body  (see  figure  1) .  Second,  the  non- 
dimensional  radii,  Sj,  S2,  connecting  the  primary  masses  to  the 
center  of  mass  are  given  by 


“  Si+S, 


where  and  S2  are  the  actual  radius  magnitudes.  The  center  of 
mass  position  Sj  as  measured  from  mj,  is  then  calculated  by 

_  _  miXO+iR^xCs.+s^)  ... 


Remembering  that 


S1+S2  =  %+m2  =  1 


then  these  parameters  may  be  redefined  by  a  single  parameter  |i. 

Si  =  i7J2  =  H  S2  =  IT?!  =  1-H  (5) 

Third,  the  non-dimensional  orbit  period  T  is  defined  as 


T  =  T 


'  G{Mj^+M2) 

GiM^+M^) 

[GCAfi+Afg) 

(51+5^2)3 

=  271 


where  G  is  the  universal  gravitation  constant.  As  a  consequence 
of  eq(6),  the  non-dimensional  angular  velocity  of  the  rotating 
coordinate  frame,  w,  may  be  derived  from  the  mean  motion  of  the 
orbiting  primaries. 


(51+^2)  3  ■ 

1 

2 

GiM^+M^) 

2 

(^1+^2)"  ■ 

G{Mj,+M2) 
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Figure  1.  Reference  Frame  for  the  Restricted  Three-Body  System 


The  Restricted  Three-Bodv  Equations  of  Motion 

In  his  treatise  on  the  problem  of  three  bodies,  Victor 
Szebehely  devotes  the  entire  first  chapter  to  the  description  and 
derivation  of  the  equations  of  motion  for  the  restricted  three- 
body  problem  (5:7-22).  The  equations  of  motion  to  be  presented 
differ  from  those  of  Szebehely,  since  his  final  result  is  a  set 
of  two  second-order  differential  equations  (3:22) .  Here,  four 
first-order  differential  equations  are  derived,  v;hich  are 
dynamically  equivalent  to  those  of  Szebehely. 

The  position  vector  for  the  third-body,  and  the  angular 
velocity  vector  for  the  rotating  coordinate  frame  can  be  directly 
observed  from  figure  1. 

r  =  +  g2'^2 

w  -  1^3  (9) 

where 

(10) 

form  a  set  of  orthogonal  unit  vectors.  The  velocity  vector  is 

-At  =  V  =  (<5-i-g2)<?i  +  (11) 

and  the  dimensionless  kinetic  energy  of  the  third-body  is  then 

r  =  =  |%[(4i-g2)"+(<5r2+gi)^]  (12) 

The  non-dimensional  potential  energy  of  is  purely 
gravitational,  and  is  given  by 


7 


GiMj^+Mz) 

.  ^2  , 

-^3  (1-H) 


+ 


^2 


(13) 


where 


^ 

^2  =  [(CTi+l-lif +(<72)"]* 

the  dimensionless  Lagrangian  is  then 


(14) 


The  conjugate  momentum  terms,  Pi  and  pj,  needed  to  complete 
the  state  vector 


x(t) 


Qi  ( t)  ■ 
Pi(t) 

<72  (t) 
P2(t) 


(16) 


are  constructed  by  differentiating  the  Lagrangian. 


Pi  =  -^  =  4-i-sri 

A  rearrangement  and  substitution  of  the  momenta  into  the 
Lagrangian  then  yields 


Thus,  the  system  Hamiltonian  is 


(18) 


H  =  -  L  =  -KPi^+PgZ)  +  Pig2  -  ^  (19) 

+  1  +2 
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a'he  equations  of  motion  for  the  restricted  third-body  may  then  be 
obtained  via  Hamilton's  equationts. 


A-  -  ^  ^ 

-  Pi  +  0-2 

~  " - 7^ . .  '  — 7^ — 

■'  '  (20) 

A  -  9/f  _  ^  ^ 

^2  Qx 

^  _  dH  _  ^  (1-J»><f2  ucr-; 

A  =  =  "Pi  " - 5 —  '  -  -r 


Periodic  Orbits  and  the  Equations  of  Variation 

Periodic  orbits  are  a  very  special  sub-set  of  restricted 
three-body  orbits.  In  general,  a  periodic  orbit  is  nothing  more 
than  an  orbit  that  closes  upon  itself  after  each  revolution,  and 
nothing  less  than  one  of  the  few  known  solutions  to  restricted 
three-body  dynamics.  Since  it  is  the  point  of  this  study  to  gain 
insight  into  nearly  periodic  orbits,  an  understanding  of  the 
periodic  case  must  first  be  obtained. 

A  periodic  orbit  will  always  return  to  it's  original  state 
after  each  integer  multiple  of  it's  period.  Because  of  this,  a 
periodic  orbit  may  be  calculated  by  iteratively  narrowing  the 
difference  between  the  initial  and  final  state  conditions. 

In  general,  periodicity  is  obtained  when 

X(0)  =  Xix)  (21) 
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In  practice,  once  a  set  of  initial  conditions  have  been 
chosen  and  the  orbit  integrated,  one  will  find  that  the  initial 
and  final  conditions  will  not  agree.  Therefore,  before  the  next 
iteration,  the  initial  conditions  must  be  adjusted  based  on  the 
error  found  in  the  final  conditions  of  the  last  integration. 

This  can  only  be  accomplished  if  dynamical  information  about 
nearby  trajectories  is  made  available.  This  information  i'-^ 
contained  w.irhin  the  equations  of  variation. 

In  vector  form,  the  equations  of  motion  may  be  ; ^written  as 


4?  =  7ix)  (22) 

Similarly,  if  we  define  the  state  of  a  nearby  trajectory  as  the 
vector  sum  of  the  current  trajectory  and  the  displacement  vector 
separating  the  two  trajectories,  then 

Xnsarby  =  X  +  ^23) 

Substituting  this  result  into  rhe  equations  of  motion,  we  get 

-§^X  *  ^(Hx)  =  7{x*6x)  (2«) 

After  expans-ion  in  a  Taylor's  series,  centered  about  the  original 
trajectory,  the  equations  of  motion  become 


_  fix  +  0|5xP 


(25) 


If  we  assume  a  first-order  expansion  and  then  subtract  the 
original  equations  of  motion,  v;e  obtain  the  equations  of 
variation. 
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In  vector  form,  Hamilton's  equations  may  be  written  as 

7(X)  =  (Z’) 


where 


0  1  0  O' 

z  =  °  0  0  (28) 

0  0  0  1 

0  0-10. 

Differentiating  Hamilton's  equations  with  respect  to  the  state 
vector  yields 

=Z-^=A(£:)  (29) 

35t  5  ap 

Where  A(t)  is  known  as  the  linearization  of  the  dynamics.  Thus, 
the  equations  of  variation  become 

(5X)  =  A(t)  5x  (30) 

expanding  A{t) , 


Ait)  = 


0 

1 

1 

0 

-^13 

-1 

0 

0 

-1 

-^33 

where 


11 


^11  = 


-3 


7*  3  7"  3 

■^1  *^2 


_  "3(gi-n)g2(l-^)  3  (<ari+l-n)  gjl* 

iil3  -  =  “  = 


=  ^13 


(32) 


^  _-3g2^a-jx)  _  3g^  ^  IzE  ^  JL 

33  _c  _q  _•»  __■» 


Since  the  equations  of  motion  form  a  set  of  first-order, 
time-varying,  linear,  differential  equations,  the  general 
solution  may  be  constructed  from  the  fundamental  set  of  solutions 
(6:61)  . 

«X(t)  =  *(t)  «X(0)  <33) 

where  ^  is  the  solution  to  the  differential  equation 


-^4(t)  =A(t)  *(t)  (34) 

at 

which  in  turn  must  be  integrated  along  with  the  equations  of 
motion. 

Upon  completion  of  an  integration,  the  solution  to  the 
equations  of  variation  may  be  constructed  as 


'5^1  (t)' 

’<^11  4>12  <1>13  4>i4 

5Pi(t) 

4>21  4>22  4>23  ^>24 

5g2<'f) 

^31  ^32  ^33  4>34 

5P2<'^) 

.4*41  4>42  4>43  4^44, 

'5gi(0) 

5pj(0) 

^QziO) 

^PziO) 


(35) 


If  in  the  selection  of  the  initial  conditions,  we  restrict  5pi(0) 
and  5g2(0)  to  zero,  then  the  orbit  must  not  only  close  upon 
itself,  but  must  do  so  intersecting  the  axis  perpendicularly. 
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The  benefits  of  this  restriction  are  three  fold.  First,  the 
number  of  initial  conditions  that  must  be  correctly  determined  is 
halved.  Second,  the  final  values  of  6pi  and  5g2  are  the  actual 
error  in  the  boundary  conditions.  Third,  only  the  portion  of 
eq(35)  that  relates  5pi('t)  and  5q-2(t)  to  5gi(0)  and  5p2(0)  is 
relevant.  Therefore,  eq(35)  may  be  simplified  to 

8Pi(x) 

and  after  matrix  inversion  to 


^21  4*24 

8gri(o) 

4*31  4*34. 

6p2(0) 

(36) 


fi<3ri(o) 

■4*ai 

4*24 

-1 

6Pi(‘c) 

8P2(0) 

.<^31 

4*34. 

(37) 


which  yields  the  necessary  relationship  between  initial  and  final 
conditions  to  find  a  periodic  orbit. 


Initial  Condition  Determination 

Since  the  equations  of  variation  are  a  linearized 

approximation  of  the  true  system  dynamics,  their  use  in  an 

iteration  scheme  is  limited.  More  to  the  point,  the  iteration 

scheme  will  only  converge  if  the  initial  conditions  chosen 

produce  a  nearly-periodic  orbit.  Thus,  a  solution  may  be 

extracted  only  if  it  has  already  been  approximated. 

The  restricted  three-body  problem  is  spanned  by  four 

dimensions,  and  is  solvable  only  by  four  exact  integrals  of 

motion.  It  was  proven  by  Henri  Poincar6,  however,  that 

The  Hamiltonian  is  the  only  analytic  integral  of  the 
motion.  If  other  so  called  'quasi  integrals'  exist. 
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they  are  not  analytic  functions  of  the  system 
coordinates,  momenta,  and  time.  (9:132) 

Therefore,  the  periodic  solution  we  seek  must  involve  the 

Hamiltonian,  and  three  of  these  'quasi  integrals' . 

Given  the  absence  of  an  analytical  approach,  it  appears  that 

periodic  regions  may  only  be  identified  through  numerical  search. 

Fortunately,  such  methods  have  been  well  developed  and  widely 

used,  given  the  availability  of  fast  and  powerful  computers. 

William  Jefferys  compiled  an  extensive  catalog  of  restricted 

three-body  phase  plots,  known  as  surface  of  section  plots  (3:1). 

These  plots  allow  the  user  to  graphically  locate  regions  of 

periodicity,  and  to  identify  sufficiently  periodic  sets  of 

initial  conditions. 

The  Surface  of  Section 

The  equations  cf  motion  as  derived  by  Jefferys  differ  from 
Szebehely's  formulation  in  two  ways.  First,  the  coordinate  frame 
is  not  centered  at  the  system  center  of  mass  point,  but  has  been 
translated  a  distance  p.  to  the  center  of  the  primary  body  of  mass 
l-[i.  (see  figure  2)  Second,  the  Jacobi  constant  C,  is  used 
instead  of  the  system  Hamiltonian,  H.  Therefore,  in  order  to 
obtain  Jefferys'  results  using  the  dynamics  already  presented, 
two  transformation  relationships  must  be  established. 

The  relationship  between  coordinates  is  simply 

(38) 

Qz=y 

while  the  relationship  between  constants  of  motion  is  a  bit  more 
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involved.  Jefferys  uses  the  Jacobi  integral  in  the  form 


where 


=  2D  -  C 


ri  =  [x^+y^]*  Tj  =  [(x+l)2+y2]* 


from  Szebehely's  equations  of  motion 

=  Pi+<3'2  ^2  =  P2-<3‘i 

so  after  substitution,  Jacobi's  integral  becomes 


(Pi+<l2)^  +  =  2 


^2 


grouping  the  terms  found  in  the  Hamiltonian  to  the  left 


H  =  ^  (p!+p|)  +  PiOTz  -  PzQi  - 
^  ^2 

=  |[-(‘3f+<22)  +  +  nr|  -  c] 


(39) 


(41) 


(42) 


(43) 


then  by  substituting  the  radius  terms 

2H  =  -(gf+gff)  +  (i-n) [ (qri-fi>^+'3f]  +  n[ (gi-|i+i)^+<3f2]  -  c  (44) 
which  simplifies  nicely  to 

2H  =  (1-jl)  -  C  (45) 

Thus,  the  transformation  between  the  derivations  of  Szebehely  and 
Jefferys  is  complete. 
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The  surface  of  section  does  not  plot  orbit  trajectories  per 

se.  Rather,  it  illustrates  the  behavior  of  all  possible 

trajectories  in  a  particular  finite  portion  of  the  phase  space. 

Because  the  Hamiltonian  is  a  constant  of  the  motion,  the 
third  body  m3  is  constrained  to  move  on  a  three  dimensional 
manifold  embedded  in  the  four  dimensional  phase  space.  If 
another  independent  integral  exists  for  this  orbit,  then  the 
third  body  would  then  be  constrained  to  move  on  a  two- 
dimensional  manifold  embedded  in  the  three-dimensional  phase 
space.  (3:6) 

The  reduction  of  dimension  from  four  to  three,  occurs  because 
specification  of  the  Hamiltonian  constant  H  and  any  three  state 
variables,  dictates  the  value  of  the  fourth.  In  this  way,  the 
state  vector  may  be  transformed  to  include  three  state  variables 
and  one  constant. 


x(t) 


<71 

Pi 

_ 

Pi 

<72 

<72 

P2. 

(46) 


The  dimension  of  the  problem  is  further  simplified  by  arbitrarily 
specifying  a  plane,  through  which  all  of  the  two-dimensional 
trajectories  of  interest  must  pass.  The  plane  chosen  by  Jefferys 
was 


XX  +  yyr  +  g^p^  -  {p^+g^)  =0  (47) 

since,  "it  represents  the  condition  for  a.',  orbit  to  have  a 
pericenter  or  apocenter  (i.e.,  the  point  nearest  to  or  farthest 
from  the  primary  l-p.)  "  (3:8).  Again,  by  algebraic  manipulation, 
the  state  vector  may  now  be  written  as 
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x(t)  = 


(48) 


Therefore,  by  specifying  a  value  for  the  Hamiltonian,  the 

arbitrary  plane,  and  the  parameter  p.,  the  dimension  of  the 

problem  is  reduced  from  four  to  two.  In  order  to  find  periodic 

regions  in  the  remaining  two-space,  Jefferys  contends 

If  a  second  integral  exists,  the  intersection  of  the  two- 
dimensional  manifold  on  which  the  particle  is  constrained  to 
move  with  the  arbitrary  plane  will  be  one-dimensional,  in 
general  (i.e.,  a  set  of  closed  curves).  On  the  other  hand, 
if  no  such  integral  exists,  then  the  intersection  will  not 
be  restricted  to  one-dimensional  sets  in  the  arbitrary 
plane.  (3:6) 

Thus,  we  now  have  a  way  to  identify  periodic  regions  in  the  phase 
space,  regions  where  two  integrals  of  motion  exist. 

Figure  3  represents  a  precessing  elliptical  trajectory  about 
primary  1-p.  By  comparison,  figure  4,  shows  which  points  from 
this  trajectory  actually  intersect  the  plane  of  interest.  In 
this  case,  the  points  of  apogee  and  perigee  constitute  the  only 
points  on  the  rotating  ellipse  that  pass  through  the  plane. 
Therefore,  the  surface  of  section  appears  as  two  circles  that 
trace  the  path  of  the  apogee  and  perigee  points,  as  the 
elliptical  orbit  processes. 

In  general,  on  most  surface  of  section  plots,  the  location 
of  the  primaries  and  mass  are  omitted.  Therefore,  one  should 
remember  that  the  plot  is  always  centered  about  the  primary  mass 
1-[L,  with  the  other  primary  mass,  [L,  located  at  the  point  (-1,0). 
The  actual  location  of  mass  is  of  little  use  and  is  never 
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Fieure  4.  The  Surface  of  Section  of  the  Elliptical  Trajectory 
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recorded.  All  that  is  really  important  here  are  the  size,  shape, 
and  location  of  closed  curves  on  the  arbitrary  plane. 

Since  it  is  the  purpose  of  these  plots  to  describe  the 
dynamics  of  a  region,  several  trajectories  must  be  integrated  and 
overlaid  to  produce  a  meaningful  surface  of  section  plot.  Figure 
5  is  a  surface  of  section  plot  for  the  sun-Jupiter  system.  The 
value  of  \L  is  small,  which  allows  %  to  orbit  the  primary  l-.U  in 
gently  perturbed  two-body  fashion.  Figure  6,  on  the  other  hand, 
is  a  good  example  of  the  highly  perturbed,  non-two-body  case. 
Here,  both  primaries  are  of  a  size  and  proximity,  that  they 
grapple  continuously  for  dynamical  control  of  On  either 

plot,  the  concentric  enclosed  island  structures  indicate  regions 
of  periodic  motion.  These  regions  are  all  centered  by  a  single 
point  that  is  a  periodic  solution. 

The  Nearlv-Periodic  Trajectory  in  Modal  Variables 

To  summarize,  the  existence  of  Jacobi's  integral  reduces  the 
dimension  of  the  phase  space  from  four  to  three.  This  is  true 
anywhere  in  the  restricted  three-body  phase  space.  Locally, 
however,  a  closed  curve  on  the  surface  of  section  plot  implies 
the  presence  of  a  second  integral.  A  second  integral  further 
confines  the  dynamics  to  a  two-dimensional  manifold,  embedded 
within  the  four-dimensional  phase  space.  The  true  dynamical 
nature  of  these  local  two-dimensional  manifolds  is  distorted  on 
the  surface  of  section  plot,  because  it  is  a  projection  of  the 
two-dimensional  manifold  onto  the  arbitrarily  defined  plane. 
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Figure  5.  The  Surface  of  Section  For  the  Sun- Jupiter  System 


Figure  6.  The  Surface  of  Section  For  a  Highly  Perturbed  System 
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Figures  7  and  8  illustrate  these  periodic  regions,  and  their 
distorted  appearance  on  a  surface  of  section  plot. 

To  construct  an  exact  representation  of  a  nearly-periodic 
trajectory  without  distortion,  a  vector  space  that  is  locally 
tangent  to  the  Hamiltonian  surface  must  be  constructed.  To  this 
end,  Wiesel  introduces  two  variable  transformations  (7:233).  The 
first,  translates  the  frame  of  reference  so  that  it  is  centered 
on  the  periodic  trajectory.  The  second,  termed  the  modal 
transformation,  orients  the  frame  with  the  local  tangent  space. 
The  resulting  transformation  maps  a  second  state  variable  into  a 
second  constant .  This  new  constant  refers  to  the  orbit  epoch 
time,  and  may  be  treated  like  the  Hamiltonian  constant  (7:236). 

These  variables  transformations  will  be  used  in  two 
different  ways.  First,  an  orbit  will  be  integrated  in  the 
original  cartesian  coordinates,  and  then  transformed  into  the 
modal  variables.  A  plot  of  this  result  will  be  used  as  an  exact 
representation  of  a  nearly-periodic  orbit.  A  wide  range  of 
displacements  will  be  tested  in  order  to  obtain  a  rough  idea  of 
the  limits  of  the  two-integral  region. 

Second,  the  original  Hamiltonian  will  be  canonically 
transformed  into  the  modal  coordinates  and  expanded  in  a  Taylor' s 
series.  Since  the  magnitudes  of  the  modal  variables  are 
typically  much  smaller  than  one  dimensionless  length,  the  higher 
order  terms  in  the  expansion  rapidly  approach  zero.  For  this 
reason,  the  Hamiltonian  expansion  will  be  truncated  after  four 
terms.  An  approximation  for  the  equations  of  motion  for  nearly- 
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Figure  7.  Magnified  Periodic  Region  of  Sun-Jupiter  SOS  Plot 
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Figure  8.  Magnified  Periodic  Region  of  Highly  Perturbed  System 

SOS  Plot 
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periodic  orbits  will  result.  The  validity  of  the  expanded 
equations  can  then  be  tested  vs.  the  exact  representation. 
The  first  transformation  begins  by  allowing 


8X  = 


figi 

fiPi 

Pi“Pio 

8^2 

Qz-Qzo 

«P2 

P2~Pz0, 

(49) 


to  constitute  a  canonical  transformation  to  a  set  of  nearly 
periodic  variables,  centered  about  a  periodic  trajectory.  Here, 
the  zero  subscript  refers  to  the  periodic  reference  trajectory, 
and  the  delta  prefix  signifies  the  new,  nearly-periodic 
variables.  A  generating  function  Fj  exists 


Fa  =  («Pi+Pio)  *  (8P2+P20)  (gi.-'^2o> 


such  that  it's  partial  derivatives  of  the  form 


Pi  = 


dQi 


^^2 

d(5p.) 


(50) 


(51) 


yield 

Pi  =  fiPi+Pio  Sg-i  =  gi-gjo  (52) 

which  reproduce  the  original  substitution  definition,  and  prove 
the  transformation  canonical. 

In  practice,  this  first  transformation  is  more  difficult 
than  it  appears.  In  order  to  subtract  the  current  nearly- 
periodic  trajectory  from  the  periodic  one,  both  trajectories  must 
exist  in  an  approximately  continuous  manner,  given  the  tiny  time 
steps  used  in  the  numerical  integration.  The  method  employed 
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here,  was  to  preserve  the  periodic  trajectory  in  a  finite  series 
of  Fourier  coefficients.  In  this  way,  barely  a  hundred  samples 
of  the  actual  trajectory  were  transformed  into  two  sets  of  fifty 
coefficients.  The  advantage  of  the  Fourier  representation  is 
that  the  coefficients  may  be  reassembled  into  the  periodic  orbit 
at  any  time  necessary.  This  is  why  the  Fourier  representation  is 
considered  continuous.  Details  of  this  digital  to  analog 
conversion  process  are  contained  in  the  Fourier  subroutine  in 
appendix  B,  and  also  in  the  book  by  Brouwer  and  Clemence  (1:109) . 

If  the  variables  describing  the  motion  of  the  third  body 
relative  to  the  periodic  orbit  are  indeed  small,  then  Wiesel 
argues  that  the  equations  of  motion  for  a  nearly-periodic 
trajectory  are  analogous  to  the  equations  of  variation  previously 
described.  In  other  words,  the  periodic  and  nearly-periodic 
trajectories  are  close  enough  that  there  is  a  linear  relationship 
between  the  two  orbits.  Thus 

-^[8x(t)]  =  A(t)  (53) 

is  a  set  of  linear,  time  periodic  differential  equations,  and 
have  the  solution 

bxit)  =9{t)  bxiO) 

According  to  Floquet,  since  A(t)  is  time  periodic,  then  the  above 
relation  may  be  decomposed  to 

5x(t)  =  F(t)  3^*^  bxiO) 

where  Af  is  a  constant  matrix,  and  F(t)  is  periodic  with  the  same 
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period  x  as  the  original  orbit.  Wiesel  then  defines  the  second 
transformation  by 

5x(t)  =  F(t)  5(t)  (56) 

where  b(t)  is  the  product  of  and  an  initial  constant. 

Because  are  the  time-varying  analogs  of  the  system 
eigenvalues,  and  F{t)  is  constructed  with  the  system 
eigenvectors,  this  transformation  is  known  as  the  modal 
transformation  (2:671).  The  new  state  vector  b(t)  is  then 
obtained  by  simple  matrix  inversion. 

5(t) 

The  orientation  of  the  unit  vectors  in  the  modal  system  is 
shown  in  figure  9.  The  matrix  differential  equation  describing 
the  periodic  changes  in  the  eigenvector  matrix  is 

-^F(t)  =A(t)F(t)  -F{t)J{t)  (58) 

dt 

where  A{t)  is  the  same  as  in  eq(31),  and  J  is  a  matrix  of  the  two 
system  Poincar^  exponents. 


JiHj  0  0 
FOj  0  0 

0  0  Re2  lai^ 

0  0  -JjT^  Re^ 


(59) 


More  information  on  the  subtleties  of  Floquet  theory  is  provided 
by  Calico  (2:672) .  The  details  proving  the  second  transformation 
canonical  are  outlined  by  Wiesel  (7:234),  where  the  reader  is 


29 


Periodic  Reference  Trajectory 


Figure  9.  Reference  Frame  in  Modal  Coordinates 
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referred  to  the  text  by  Pars  (4:453-483)/  for  a  full  explanation 
and  proof. 


The  Restricted  Three-Body  Perturbation  Solution 

After  canonical  transformation  to  the  modal  coordinates,  the 
new  Hamiltonian  minus  a  pure  function  of  time  is 

K(B)  =  H(B)  (60) 

since 


dt 

doesn't  contribute  to  the  equations  of  motion. 

Expanding  the  new  Hamiltonian  in  a  Taylor' s  series  produces 


4  4 


kcb)  =  H(o)  +  y  y  zME.  B  B  + . . . 

where  b=0  centers  the  expansion  about  the  periodic  trajectory. 
Alternatively,  this  expansion  may  be  more  compactly  written  using 
tensor  notation. 


kCB)  =  H(0)  +  +  •  •  • 


The  first  term  in  the  expansion  is  the  Hamiltonian  for  a 
periodic  orbit,  and  is  a  constant.  The  second,  or  linear  term  is 
identically  zero,  because  it  describes  the  motion  of  the  periodic 
trajectory  with  respect  to  itself.  The  third,  or  quadratic  term 
is  the  Floquet  problem,  and  becomes  a  constant  coefficient, 
linear  system  in  the  new  variables.  Since  the  magnitude  of  the 
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modal  state  vector  is  very  small  compared  to  one,  the  expansion 
is  truncated  after  the  fourth  term. 

The  first  two  elements  of  the  modal  state  vector,  bi  and  bz, 
are  the  only  two  variables  in  the  modal  space.  Elements  bj  and 
jb^,  are  both  constants,  and  represent  a  change  in  the  orbit  epoch 
(bj) ,  and  a  change  in  the  Hamiltonian  surface  (b^)  .  Since  we  are 
free  to  arbitrarily  choose  both  of  these  values,  they  are  both 
set  to  zero.  The  dimension  of  the  new  Hamiltonian  is  then 
reduced  from  four  to  two. 

KS)  =  mo)  ^  I  S'-  “  °]  5  *  -i  5,  3,  3,  («<) 

where  CD  is  the  imaginary  portion  of  the  non-zero  Poincar6 
exponent,  and  the  modal  state  vector  is 

5  = 

Expanding  the  third-order  tensor,  the  new  Hamiltonian  becomes 

K(B)  =  H(0)  +  (bi+bf)  (0  +  jbiOi  -  (66) 

where 


32 


«!  ( t)  “  — 

1 

“  ~^[^ijk^j.i^2j^k  ■*■  ^ijk^i^lj^k  '*'  ^Ijk^zi^j^k] 
~  ^^ijk^l^2j^2k 


which  are  periodic  functions  of  time  alone.  and  Ai,  are 

functions  of  the  periodic  trajectory  and  are  independent  of  the 
modal  variables.  A.,  represents  a  particular  eigenvector  element 
from  the  eigenvector  matrix.  Here/  the  index  order  is  reversed 
such  that  i  represents  a  particular  eigenvector,  and  ^  indicates 
which  element  from  that  eigenvector  is  needed.  Using  Hamilton's 
equations,  the  equations  of  motion  for  nearly  periodic  orbits  are 


=  ojbj  -  +  2iJib2a3  -  3bla^ 

-^i?2  =  -(oiJi  -  jbfaj  +  2biZ72«2  “ 


(68) 
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IV.  Hardware  -  Software 


This  study  is  a  collection  of  several  programs,  coded  in 
standard  Fortran  77,  and  executed  on  an  ELXSI  6420  Super-Mini 
computer.  All  of  the  software  is  contained  within  appendix  B, 
and  is  made  up  of  six  different  programs  and  their  associated 
subroutines.  The  final  computer  outputs  are  a  collection  of  data 
files  composed  of  ordinate  and  abscissa  values.  The  figures 
presented  below  were  created  by  simply  plotting  the  unformatted 
data  values,  (x,y). 
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V.  Numerical  Technique 


Searching  the  Phase  Space 

Before  a  periodic  trajectory  can  be  extracted,  periodic 
regions  in  the  phase  space  must  first  be  identified  using  program 
SECTION.  Isolating  a  periodic  region  and  a  useful  set  of  initial 
conditions  requires  a  cervain  amount  of  trial  and  error  without 
prior  knowledge  of  the  phase  space. 

Initial  Conditions 

Once  a  periodic  region  has  been  located  on  a  surface  of 
section  plot,  a  set  of  initial  conditions,  (xO,yO),  a  value  for 
the  Hamiltonian  constant,  and  a  value  for  the  parameter  |l  can  be 
determined.  The  orbital  period  must  also  be  estimated  by 
observing  the  number  of  integration  steps  needed  for  the  orbit  to 
return  to  it's  original  state,  and  multiplying  this  value  by  the 
time  step  (time/step) . 

Extracting  a.  Periodic  Trajectory 

Using  the  initial  conditions  found  above,  a  periodic 
trajectory  is  extracted  from  the  periodic  region  by  the  program 
PERIOD.  The  initial  value  of  the  orbital  period  must  be  adjusted 
in  order  to  find  the  period  that  corresponds  to  the  Hamiltonian 
surface  of  interest.  There  is  an  inverse  relationship  between 
the  orbit  period  and  the  Hamiltonian  surface  constant.  Upon 
convergence,  the  program  will  calculate  the  eigenvectors, 
eigenvalues,  and  Poincar6  exponents  of  the  linearized  system. 
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Storing  the  Periodic  Trajectory 

In  order  to  preserve  the  periodic  trajectory  in  nearly 
continuous  form/  the  state  vector  and  eigenvector  matrix 
previously  calculated  must  be  fed  into  program  FLOQUET/FOURIER. 
This  program  numerically  integrates  the  periodic  trajectory/  and 
the  state  transition  matrix.  Routine  samples  of  the  state  vector 
and  the  eigenvectors  of  the  state  transition  matrix  are  taken  and 
placed  in  temporary  storage-  Upon  completion  of  the  integration, 
the  stored  values  are  fed  into  a  Fourier  conversion  subroutine. 
One  hundred  of  the  possible  thousand  integration  values  are 
converted  into  two  sets  of  fifty  Fourier  coefficients. 

Storing  the  Periodic  Hamiltonian  Coefficients 

The  program  FLOQUET/HAMILTONIAN  requires  the  same  input  as 
program  FLOQUET/FOURIER.  Here/  the  periodic  trajectory  is 
integrated  and  sampled  as  before.  The  periodic  Hamiltonian 
coefficients  are  then  calculated  by  a  complex  series  of  vector 
multiplications.  As  before/  the  values  are  temporarily  stored 
until  the  integration  is  complete.  Another  call  to  the  Fourier 
subroutine  produces  the  desired  set  of  coefficients. 

The  Exact  Nearly-Per iodic  Trajectory 

The  exact  nearly-periodic  trajectories  are  created  by  the 
program  EXACT.  The  orbit  is  integrated  in  the  original 
variables/  before  the  state  vector  is  routinely  extracted  and 
transformed  into  the  modal  variables.  Three  separate  data  files 
result.  One  each  for  the  two  constants  created  by  the  modal 
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transformation  vs.  time/  and  one  for  the  remaining  two  variables 
plotted  together.  The  only  inputs  required  include  the  Fourier 
coefficients  of  the  periodic  orbit/  and  a  value  for  the  initial 
displacement  off  the  periodic  center.  The  output  is  an 
unformatted  data  file/  whose  (X/y)  entries  may  be  plotted 
directly. 

The  Expanded  Nearlv-Periodic  Trajectory 

Here/  the  new  set  of  equations  derived  from  the  expanded 
Hamiltonian  are  integrated.  The  inputs  required  include  the 
Fourier  coefficients  that  represent  the  periodic  Hamiltonian 
coefficients/  and  the  same  initial  displacement  used  in  program 
EXACT  converted  to  modal  coordinates.  Program  EXACT  provides 
these  values.  The  output  is  a  data  file  containing  both  elements 
of  the  state  vector/  sampled  during  the  integration. 

Numerical  Analysis 

There  are  two  distinct  observations  to  made  in  this 
analysis.  First/  several  exact  trajectories  will  integrated  and 
grouped  according  to  their  initial  displacements  from  the 
periodic  trajectory.  A  rough  idea  of  the  behavior  of  nearly- 
periodic  orbits  vs  initial  displacement  can  be  obtained.  Also/ 
the  validity  of  the  two  integral  assumption  will  be  monitored  by 
plotting  both  constants  of  motion  vs  time. 

The  second  set  of  observations  will  be  made  by  overlaying 
several  expanded  orbits  on  top  of  the  exact  ones.  This  step  is 
critical  in  determining  the  proper  truncation  limit  for  the 
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expanded  Hamiltonian.  Both  sets  of  observations  will  be 
accomplished  twice.  Once  for  the  Sun- Jupiter  system,  and  once 
for  the  highly  perturbed  system. 
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VI.  Results  and  Discussion 


Nearlv-Periodic  Orbits  in  Modal  Variables 

Figure  5  is  the  surface  of  section  plot  for  the  Sun-Jupiter 

system.  Here,  |i  is  very  small,  and  as  one  might  expect,  this 

system  could  well  be  described  using  classical  two-body 

perturbation  theory,  were  it  not  in  the  vicinity  of  a  resonance. 

Proximity  to  a  resonance,  means  that  the  same  relative 
configurations  repeat  in  the  same  order.  This  gives  an 
otherwise  small  gravitational  force  a  chance  to  produce  a 
relatively  large  effect.  (9:98) 

By  comparison,  figure  6  represents  a  highly  perturbed  dynamical 
system,.  These  are  only  two  of  an  infinite  number  of  choices,  but 
were  chosen  because  of  the  plainly  visible  closed  curves  present 
on  both.  Figures  7  and  8,  magnify  these  regions  where  two 
integrals  of  motion  are  present. 

To  this  point,  all  of  the  figures  presented  have  been  in  the 
coordinate  frame  presented  in  figure  2.  From  this  point  on,  all 
of  the  figures  will  be  referenced  to  the  modal  unit  vectors 
described  in  figure  9,  and/or  a  time  axis  measured  in  orbital 
radians. 

The  exploration  of  nearly-periodic  three-body  trajectories, 
begins  with  figure  10.  Here,  orbits  very  close  to  the  Sun  in  the 
Jupiter-Sun  system  are  plotted.  As  was  expected,  a  set  of 
concentric  circles  is  present.  The  initial  displacement  off  the 
periodic  trajectory  for  each  orbit  is  recorded  as  a  fraction  of 
ji.  This  plot  depicts  the  limiting  region  where  two  integrals  are 
said  to  exist.  Any  distortion  of  the  trajectories  indicates  the 
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Figui*e  10.  Neariy-Periodic  Orbits  in  the  Vicinity  of  Two  Exact 
Integrals  of  Motion  (Sun-Jupiter  System) 
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gradual  dissolution  of  the  second  constant.  Figure  11  shows  the 
behavior  of  the  two  integrals  for  the  outer  trajectory  of  figure 
10,  plotted  on  the  same  scale.  The  small  periodic  displacements 
on  figure  11  correspond  to  the  distortions  on  figure  10.  Since 
straight  lines  and  perfectly  round  circles  are  all  that  exist  for 
the  trajectories  inward  from  this  point,  they  were  considered 
uninteresting  and  left  alone. 

The  marginal  case,  where  the  second  integral  of  motion  can 
no  longer  be  assumed  constant,  is  pictured  in  figure  12.  Here, 
the  initial  displacements  off  the  periodic  trajectory 
are  increased  by  an  order  of  magnitude.  The  circular  orbits 
appear  to  decompose  into  five  separate  trajectories. 

Surprisingly,  however,  each  orbit  contir.w-.es  to  close  upon  itself. 
Figure  13  shows  the  marginally  constant  nature  of  the  second 
integral. 

The  extreme  case,  where  Jacobi's  integral  is  the  only 
constant  in  the  system,  is  plotted  in  figure  14.  Unbelievably, 
these  orbits  about  the  sun  continue  to  close.  This  is  very 
interesting  since  the  apogee  of  the  orbit  extends  all  the  way  out 
to,  and  even  past  Jupiter!  Figure  15  plots  the  third  state 
variable,  formerly  the  second  integral,  and  Jacobi's  constant. 

The  next  series  of  plots,  figures  16-21,  represent  similar 
cases  to  the  ones  presented,  only  now  for  the  system  where  M-  is 
equal  to  a  third.  Once  again,  the  first  plot  in  the  series, 
figure  16,  illustrates  the  nicely  concentric  nature  of  the 
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Figure  11.  Epoch  and  Hamiltonian  Constants  vs.  Time  in  the 
Vicinity  of  Two  Exact  Integrals  of  Motion  (Sun-Jupiter  System) 
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Figure  12.  Neariy-Periodic  Orbits  in  a  Transition  From  Two  Exact 
Integrals  of  Motion  to  One  (Sun-Jupiter  System! 


43 


Figure  13.  Epoch  and  Hamiltonian  Constants  vs.  Time  in  a 
Transition  From  Two  Exact  Integrals  to  One  ( Sun-Jupiter  System) 


44 


b4 


I  j  /i  =  XJ00953SS 

i  I 

I  i  !  _  7  -1  A 

I  I  u  ~  u  I  \  \J 


I  I  -  , 

I  I  >u  =  All 

i  I 


Figure  15.  Epoch  and  Hamiltonian  Constants  vs.  Time  in  the 
Absence  of  a  Second  Integral  of  Motion  ( sun- Jupiter  System) 
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Figure  16.  Neariy-Periodic  Orbits  in  the  Vicinity  of  Two  Exact 
Integrals  of  Motion  (Highly  Perturbed  System) 
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nearly-periodic  orbit  in  the  presence  of  both  integrals. 

Figure  17  confirms  that  there  are  indeed  two  constants  of  motion 
present. 

In  figure  18,  the  phenomenon  of  chaotic  motion  is  first 
seen.  As  we  step  further  away  from  periodicity,  the  orbits 
develop  increasing  numbers  of  odd  twists  and  turns.  The  outer 
three  orbits  fail  to  close  at  all.  Figure  19  is  very 
interesting,  since  for  the  first  time,  both  integrals  of  motion 
appear  to  be  moving.  Since  Jacobi's  constant  has  been  proven  to 
exist  everywhere  in  the  phase  space,  then  figure  19  indicates 
that  the  tangent  space  is  no  longer  aligned  with  the  Hamiltonian 
surface . 

Figure  20  is  an  excellent  example  of  a  formerly  well  behaved 
dynamical  system  marching  off  to  chaos.  The  inner  two 
trajectories  are  still  recognizable  as  orbits  that  nearly  close. 
The  outer  two  trajectories  no  longer  describe  an  orbit.  The 
onset  of  chaotic  motion  is  well  presented  in  figure  21.  For 
approximately  half  an  orbit,  the  second  integral  remains  quasi 
constant,  and  the  coordinate  frame  anchored  on  the  Hamiltonian 
surface.  As  time  progresses,  the  epoch  constant  becomes  the 
epoch  variable,  and  the  relationship  between  coordinate  reference 
and  the  Hamiltonian  surface  is  destroyed.  Total  chaos  ensues. 
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Figure  17.  Epoch  and  Hamiltonian  Constants  vs.  Time  in  the 
Vicinity  of  Two  Exact  Integrals  of  Motion  (Highly  Perturbed 

System ) 
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Figure  18.  Nearly-Periodic  Orbits  in  a  Transition  From  Two  Exact 
Integrals  of  Motion  to  One  (Highly  Perturbed  System) 
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Figure  19.  Epoch  and  Hamiltonian  Constants  vs.  Time  in  a 
Transition  From  Two  Exact  Integrals  to  One  (Highly  Perturbed 

System ) 
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Figure 


20.  Neariy-Periodic  Orbits  in  the  Absence  of  a  Second 
Integral  of  Motion  (Highly  Perturbed  System) 
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Figure  21.  Epoch  and  Hamiltonian  Constants  vs.  Time  in  the 
Absence  of  a  Second  Integral  of  Motion  (Highly  Perturbed  System) 
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The  Expanded  Approximation  Vs.  The  Exact  Case 

To  best  illustrate  the  accuracy  of  the  approximated 
equations  of  motion,  they  have  been  overlaid  with  a  plot  of  the 
exact  case.  Figure  22  compares  two  of  the  trajectories  from 
figure  10,  with  the  same  two  trajectories  from  the  expanded  case 
The  inner  orbits  overlay  in  round  concentric  fashion  as  expected 
The  outer  overlay  represents  the  limit  of  the  highly  correlated 
region.  In  figure  23,  the  trajectories  are  very  marginally 
agreeable,  while  in  figure  24,  the  trajectories  only  meet  at 
nodal  points. 

For  |i  equal  to  a  third,  similar  behaviors  may  be  observed. 
Figure  25  represents  the  highly  correlated  region.  Figure  26, 
however,  shows  a  large  difference  between  the  expanded  and  exact 
cases.  This  case  is  still  marginal,  since  both  represent  nearly 
closed  orbits  of  the  same  size.  In  figures  27,  the  expanded 
trajectories  don't  show  any  indication  of  the  chaotic  motion 
present  in  the  exact  case. 
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Figure  22.  Overlay  of  Exact  and  Expanded  Nearly-periodic 
Trajectories  in  the  Vicinity  of  Two  Integrals  of  Motion  (Sun 

Jupiter  System) 
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Figure  23.  Overlay  of  Exact  and  Expanded  Nearly-periodic 
Trajectories  in  a  Transition  From  Two  Exact  Integrals  of  Motion 

to  One  (Sun- Jupiter  System) 


Figure  24.  Overlay  of  Exact  and  Expanded  Neariy-periodic 
Trajectories  in  ther  Absence  of  a  Second  Integral  of  Motion  ( Sun- 

Jupiter  System) 
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Figure  25.  Overlay  of  Exact  and  Expanded  Nearly-periodic 
Trajectories  in  the  Vicinity  of  Two  Integrals  of  Motion  (Highly 

Perturbed  System) 
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Figure  26.  Overlay  of  Exact  and  Expanded  Nearly-periodic 
Trajectories  in  a  Transition  From  Two  Exact  Integrals  of  Motion 

to  One  (Highly  Perturbed  System) 
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ure  27.  Overlay  of  Exact  and  Expanded  Nearly-periodic 
Trajectories  in  the  Absence  of  a  Second  Integral  of  Motion 

(Highly  Perturbed  System) 
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VII.  Conclusions  and  Recommendations 

The  Modal  Transformation  and  the  Limits  of  the  Tangent  Space 

The  transformation  into  modal  variables  does  indeed  map  the 
state  vector  into  a  set  of  two  variables  and  two  constants  in  a 
small  region  surrounding  the  periodic  orbit.  In  the  Sun-Jupiter 
system,  the  radius  of  displacement  separating  the  regions  where 
the  epoch  was  or  wasn't  constant,  appears  to  be  about  10%  of  |i  or 
approximately  l.OE-4.  The  Hamiltonian  constant  remained  so,  for 
every  displacement  attempted.  This  is  interpreted  to  mean  that 
the  Hamiltonian  surface  is  relatively  flat  over  a  large  area, 
which  permits  an  extended  region  of  alignment  between  the 
Hamiltonian  surface  and  the  tangent  space.  Given  the  emergence 
of  epoch  time  as  a  variable,  this  technique  could  be  applied  to 
two-body  systems  near  a  resonance. 

In  the  highly  perturbed  case,  the  apparent  two-integral 
region  occurred  within  1%  of  p..  While  this  is  a  much  smaller 
percentage  of  p.  than  in  the  previous  case,  the  actual  size  of  the 
region  is  3.0E-3.  In  terms  of  a  ratio,  this  region  is  35  times 
as  large  as  the  region  in  the  Sun-Jupiter  system.  Although  this 
result  was  not  expected,  it  makes  sense  in  terms  of  the  location 
of  the  center  of  mass  of  the  three-body  system.  In  the  Sun- 
Jupiter  system,  the  mass  center  is  nearly  on  the  sun,  where  as 
the  mass  center  for  the  highly  perturbed  case  is  a  third  of  the 
distance  separating  the  primaries  away  from  the  primary  l-^i. 

As  before,  a  transition  of  the  epoch  constant  to  the  epoch 
variable  was  observed  with  larger  displacements  from  the  periodic 
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trajectory.  In  this  case,  however,  the  Hamiltonian  constant 
began  to  move  as  well.  Both  observations  coincided  with  the 
transition  to  chaos.  It  appears  that  the  Hamiltonian  surface  and 
the  tangent  space  skew  much  more  abruptly  than  before.  Unlike 
the  Sun-Jupiter  system,  the  useable  portion  of  the  phase  space  is 
strictly  limited  to  the  region  where  both  integrals  exist. 


The  Expanded  Vs .  the  Exact  Equations  of  Motion 

The  exact  and  expanded  versions  of  the  nearly-periodic 
trajectories  correlate  completely  in  the  presence  of  both 
constants  of  motion.  For  the  Sun-Jupiter  case,  within  10%  of  |i, 
and  for  the  highly  perturbed  case,  within  1%  of  |l.  These  limits 
were  determined  by  simple  observation  of  the  figures  presented, 
and  are  not  intended  to  be  exact.  In  these  cases,  however,  the 
extra  terms  in  the  expanded  Hamiltonian  equations  of  motion 
eq(68),  other  than  the  first,  are  unnecessary.  These  extra  terms 
are  either  zero,  or  very  small,  and  fail  to  accurately  model  the 
transition  of  the  epoch  variable.  Therefore,  if  the  desired 
trajectories  are  strictly  limited  to  the  two-integral  regions, 
then  the  equations  of  motion  could  be  further  reduced  to 


4  =  ^^2 
4  = 


(69) 


which  are  the  equations  of  a  harmonic  oscillator.  Since  this 
system  can  be  solved  in  closed-form,  no  integration  would  be 
necessary.  The  programs  EXPANDED  and  FLOQUET/HAMILTONIAN  could 
be  eliminated  completely. 
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Conversely,  if  the  nearly-periodic  trajectories  outside  of 
the  two-integral  region  are  to  be  explored,  then  more  terms  must 
be  maintained  in  the  Taylor's  expansion  of  the  new  Hamiltonian. 
Eventually,  by  comparison  of  the  new  expanded  trajectories 
against  the  exact  ones,  the  proper  truncation  limit  may  be 
determined.  Once  this  has  been  accomplished  the  perturbation 
solution  may  then  be  constructed. 

Either  way,  a  solution  for  nearly-periodic  orbits  can  be 
derived  in  closed-form.  The  entire  system  would  be  modelled  as  a 
harmonic  oscillator  with  a  forcing  term,  and  no  further 
integration  would  be  necessary.  The  final  task  would  then  be  to 
derive  a  functional  relationship  between  the  system  parameters 
and  integrals  of  motion,  and  the  maximum  allowable  initial 
displacement  from  the  periodic  orbit. 

{b,{0)  =  f(b„b„\i) 
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Appendix  A:  Code  Validation  and  Error  Determination 

To  ensure  that  the  computer  programs  used  in  this  study  were 
correct,  a  large  amount  of  testing  and  cross  checking  was 
required  during  code  development.  Wherever  possible,  hand  checks 
were  performed.  Unfortunately,  given  the  numerical  complexity  of 
most  of  these  programs,  few  were  actually  validated  that  way. 

The  Surface  of  Section 

It  was  very  easy  to  determine  the  accuracy  of  this  program. 
Assuming  that  the  catalog  of  SOS  plots  published  by  Jefferys  was 
correct,  the  program  was  complete  when  the  pictures  matched. 

This  step  validated  the  equations  of  motion  for  the  restricted 
three-body  problem,  as  well  as  the  surfacing  technique. 

Periodic  Orbit  Determination  Program 

The  simple  beauty  of  a  periodic  orbit  is  that  it  returns  to 
it's  original  conditions  after  each  orbit.  The  state  transition 
matrix  was  checked  by  taking  numerical  derivatives,  and  it's 
eigenvectors  using  a  linear  algebra  software  package.  The  last, 
best  test  was  a  hand  check  of  the  converged  state  values 
substituted  into  the  equations  of  motion. 

Floquet /Fourier  and  Floquet /Hamiltonian 

Both  of  these  programs  integrated  the  equations  of  motion 
already  validated,  using  the  periodic  initial  conditions,  also 
already  checked.  Even  though  the  individual  outputs  were 
different,  the  methods  used  were  the  same.  Periodic  state  or 
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Hamiltonian  information  was  extracted  at  fixed  intervals  during 
the  integration.  The  only  check  available  here  was  to  scrutinize 
the  general  periodic  trends  in  the  data.  Upon  completion  of  the 
integration,  the  discrete  periodic  information  was  converted  into 
Fourier  coefficients.  These  in  turn  were  checked  by  simply 
reconstructing  the  known  trajectory  from  the  coefficients. 

The  Program  Exact 

Here,  the  only  process  yet  to  be  checked  was  the  matrix 
inversion  from  the  intermediate  set  of  coordinates  to  the  modal 
ones.  This  was  accomplished  via  the  linear  algebra  package.  The 
two  constants  were  checked  for  any  periodic  modulation  that  may 
have  been  caused  by  insufficient  sampling  of  the  periodic 
trajectory.  In  these  cases  the  number  of  Fourier  samples  was 
simply  increased. 

The  Program  Expanded 

The  only  test  left  to  perform,  was  to  measure  the  size  and 
shape  of  the  output  compared  to  the  exact  case,  using  identical 
inputs . 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  PROGRAM  SECTION 

c 

c  PURPOSE:  Creates  a  surface  of  section  plot  file 

c 

c  SUBROUTINES:  HAMING.F 

c  RHSl.F 

c  H.F 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

program  sect 

problem  commons 

common  /data/  xmu,xmua 
common  /lam/  xlambda(4) 

common  /ham/  t,x(20, 4) ,  f  (20, 4)  ,err  (20)  ,nn,hh,mode 

variable  declarations 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
character*10  filnaml, filnam2 

dimension  xlambda(4)  ,x(20, 4) ,  f  (20, 4) , err  (20) ,  rdotv  (4) 

input  data 

read(*,*)  xmu,xmua 
read(*,*)  hh,tmax 
npts  =  dint (tmax/hh) 
read(*,*)  xnot,ynot 
read(*,*)  xjac,  syn 
read(*,*)  filnaml 
read{*,*)  filnam2 

open  output  files 

—  filnaml  is  the  plotfile 

—  filnam2  is  a  general  output  file 

open  (2,  fileofilnaml,  status»'un)cnown' ) 
open (3, file=filnam2, status^' unknown' ) 

write (3,*)  'rau=',xmu,'  l-mu=',xmua 
write (3, *j  'x0=' ,xnot, 'y0=' ,ynot 
write(3,*)  ' jacobi  const=',xjac,'timestep=',hh 
write (3, •) 

write (3,*)  'npts  value  at  n*period' 

mode  =  0 
nn  =  4 
nxt  =  0 
t  =  O.dO 

get  ql,pl,q2,p2  for  given  xO,yO,  and  jacobian 

ql  •>  xnot  +  xmu 
q2  =  ynot 

xham  =  (xmu*xmua-xjac) /2.d0 
rl  =  ((ql-xmu)»*2.d0  +  q2**2.d0) ‘‘.SdO 
r2  “  ( (ql+xmua) **2.d0  +  q2**2.d0) **.5d0 
d  =  xham  +  xmua/rl  +  xmu/r2 
g  =  q2/ (ql-xmu) 
a  =  g*g  +  l.dO 

b  =  -2.d0* (g*g»xmu  +  g*g2  +  gl) 
c  “  (g»xmu)  **2.d0  +  2.d()*g*q2*xmu  -  2.d0»d 
p2  »  (-b+syn*(b*b-4.d0»a*c)**.5d0)/(2.d0*a) 
pi  =  g*(xrau-p2) 

initial  conditions 

x(l,l)  =  ql 
x(2,l)  =  pi 
x(3,l)  =  q2 
x(4,l)  =  p2 

initialize  haming 


OOO  0000000  oooooo  oooooo  ooo 


call  hamlng{nxt) 

turn  off  second  EOM  eval 

nxt  -nxt 

lf(nxt  .ne.  0)  go  to  499 
stop  99 

499  continue 

integration  loop 

do  500  i  «  l,npts 

permute  indices 

nm3  =  nm2 
nm2  =  nml 
nml  =  nxt 

integrate  orbit,  haming  permutes  nxt 

call  haming(nxt) 

calculate  r  dot  v 

qld  «  f(l,nxt) 
q2d  =•  f(3,nxtj 

rdotv(nxtj  =  (x  (1,  nxt) -xmu)  *qld  +  x(3,nxt)*q2d 

check  for  peri/apoapse  crossing 

if  (tdotv (nxt) *rdotv{nml) .gt.O.dO)  go  to  500 

crossing  has  occured! ! ! 
interpolate  to  crossing  time 

frac  =  -rdotv(nxt) /{  rdotv(nxt)  -  rdotv{nml)  ) 
qlc  =  -frac*x(l,nml)  +  (l.dO  +  frac)  *x(l,  nxt) 
q2c  =  -frac»x (3,nml)  +  (l.dO  +  frac) ‘x (3, nxt) 
xcross  »  qlc  -  xmu 
ycross  «  q2c 

compute  conjugate  momenta  pic  and  p2c  for  qlc  and  q2c 

rlc  =  { (qlc-xmu)**2.d0  +  q2c»*2.d0) **.5d0 
r2c  =■  ( (qlc+xmua)»*2.d0  +  q2c*»2.d0)  *».5d0 
dd  =  xham  +  xmua/rlc  +  xmu/r2c 
gg  =  q2c/ (qlc-xmu) 
aa  =  gg*gg  +  l.dO 

bb  =  -2.d(5*  (gg*gg»xrau  +  gg»q2c  +  qlc) 
cc  =  (gg*xmu) »*2.d0  +  2.d0*gg*q2c»xmu  -  2.d0*dd 
p2c  =  (“bb+syn* {bb*bb-4.d0*aa*cc) **.5d0) /  (2.d0*aa) 
pic  ”  gg* (xmu-p2c) 

write{3,*)  l,x(l,nxt) 
write  (2,»)  xcross, ycross 

500  continue 

close (2) 
close (3) 
stop 
end 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  PROGRAM  PERIOD 

c 

c  PURPOSE:  Calculate  periodic  trajectory  via  iterative 

c  integration  technique.  Once  found,  determines 

c  eigenvalues,  eigenvectors,  and  Poincare  exponents 

c  for  periodic  trajectory  at  t°0. 

c 

c  SUBROUTINES:  HAMING.F 

c  RHSl.F 

c  H.F 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


program  po 

common  /lam/  xlambda(4) 
common  /data/  xmu,xmua 

common  /ham/  t,  x(20,4),f(20,4),err(20),nn,hh,mode 
local  variables 


implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

dimension  xlambda(4)  ,x(20,4) ,  f  (20, 4) , err  (20)  ,cerr  (2, 1)  ,b(2,2) 

dimension  phi  (4,4)  ,xxx(l0)  ,w)c(50)  ,xww(2) ,  rvec(2,16) 

dimension  tvec(16,2)  ,alpha(4)  ,tnvec(16,2)  ,xreal  (4)  ,ximag(4) 

complex*16  w  (4) ,  vec  (4, 4)  ,ww 

equivalence  (ww,xww) 

equivalence  (vec, rvec) 

character*10  filnaml, fllnam2 

read  input  data 

read(*,*)  xmu,xmua 
read(*,*)  period, npts 
read(*,*)  x0,y0 
read(*,*)  tol,maxit 
read(*,»)  xjacob 
read(*,*)  filnaml 
read(*,*)  filnam2 

calculate  tlmestep 
hh  a  period/ (dble (npts) ) 

calculate  ql,pl,q2,p2  for  given  x0,y0,  and  jacobian 

ql  =  xO  +  xmu 
q2  yO 

xham  "  (xrau*xmua-xjacob) /2.d0 

rl  •»  ( (ql -xmu)  **2.30  +  q2**2.d0)  “.SdO 

r2  =  ( (ql+xmua) **2.d0  +  q2**2.d0) •♦.5d0 

dd  a  xham  +  xmua/rl  +  xmu/r2 

gg  =  q2/ (ql-xmu) 

aa  =  gg*gg  +  l.do 

bb  =  -2.d0* (gg*gg*xmu  +  gg*q2  +  ql) 
cc  =  (gg*xmu) **2.d0  +  2.d0*gg*q2*xmu  -  2.d0*dd 
p2  =  (-bb+(bi5*bb-4.dO*aa*cc)**.5dO)/(2.dO*aa) 
pi  »  gg*(xmu-p2) 

echo  inputs  to  output  file 

open  (3,  file=filnaml,  status='un]cnown' ) 

write (3,*)  'xmu  =  ',xmu, '  xraua  =  ',xmua 
write (3,*)  'orbit  period,  npts  ', period, npts 
write (3,*)  'timestep  ',hh 
write(3,*)  'qlO  =«  ',ql,'  plO  =  ',pl 
write(3,*)  'q20  =  ',q2,'  p20  =  ',p2 
write (3,*)  'tol,  raaxit',toi,maxit 

begin  iteration  loop 
do  1000  iter  °  l,maxit 


set  up  initial  state 


x(l,l) 

x(2,l) 

x(3,l) 

x(4,l) 


P2 


write  progress 
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000  000  000  000000  000  000  000  000  000  ooo 


c 


101 

100 


write  (3,*) 

write (3,*)  ' iteration' ,  iter 
write (3,*)  'ql»',ql,'  pi-', pi 
write{3,»)  'q2-',q2,'  p2=',p2 

initialize  phi  matrix 

do  100  i  -  1,4 
do  101  j  -  1,4 
ij  =  4*i+j 
x(ij,l)  =  O.dO 
continue 
x{5*i,l)  -  l.dO 
continue 


initialize  integration  constants 

mode  =  1 
nn  —  20 
nxt  =  0 
t  =  O.dO 


initialize  haming 

call  haming (nxt) 

if  (nxt  .ne.  0)  go  to  499 

write  (*,»>  'failure  to  initialize' 
stop  99 
499  continue 


integration  loop 

do  500  i  =  l,npts 
call  haming  (nxt) 
500  continue 


extract  error  vector 

cerr(l,l)  =  -x(2,nxt) 
cerr(2,l)  =  -x(3,nxt) 

extract  correction  matrix 


b(l,l) 

b(l,2) 

b(2,l) 

b(2,2) 


°  x(9,nxt) 

-  x(12,nxt) 

-  x(13,nxt) 

-  x(16,nxt) 


calculate  state  corrections 


call  leqt2f (b, 1, 2, 2, cerr, idig, xxx, ier) 
add  in  corrections 


ql  =  ql  +  cerr  (1,1) 
p2  »  p2  +  cerr (2,1) 

check  for  convergence 

lend  =  0 

if  (dabs  (cerr  (1,1) )  .gt.  tol)  lend  =  1 
if  (dabs  (cerr  (2,1) )  .gt.  tol)  lend  =  1 
if  (lend  .eq.  0)  go  to  2000 

continue 


maximum  iterations  exceeded  without  convergence 

write  (*,*)  'Iteration  Limit  Exceeded' 
stop 

continue 

converged  processing 
write  (3, ») 

write (3,*)  'program  converged  in', iter, 'iterations' 
write  (3,*) 

write (3,*)  'converged  state  values' 
write (3,*)  'ql-' ,x(l,nxt) ,  'pl»'  ,x(2,nxt) 
write (3, •)  'q2-',x(3,nxt),  'p2»',x(4,nxt) 
writc(3,  •) 

write (3,*)  'surface  of  section  coordinates' 
writc(3,*)  'X  »',x(l,nxt)-xmu, 'y  =',x{3,nxt) 
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000  OOOOOO  OOO  OOO  O  000  ooooo 


write  (3,*) 
c 

c  compute  hamiltonian/ jacobian 

ql  »  x(l,nxt) 
pi  =  x{2,nxt) 
q2  “  x(3,nxt) 
p2  =  x(4,nxt) 

rl  =  dsqrt ( (ql-xmu) **2.dO  +  q2**2.d0) 
r2  =  dsqrt { (ql+xmua) **2. do  +  q2**2.d0) 

xham  =  .SdO* {pl*pl+p2*p2)  +  pl»q2  -  p2»ql  -  xmua/rl  -  xmu/r2 
xjac  ”  xrau*xmua  -  2.d0»xham 

write(3,*)  'ham  •=' , xham, ' jac  =',xjac 

extract  phi 

do  2005  i  =  1,4 
do  2005  j  =  1,4 

phi(i,j)  =  x(4»i+j,nxt) 

005  continue 

compute  eigen  values  and  vectors  of  phi 

call  eigr f (phi , 4 , 4 , 2, w, vec, 4 , wk, ier ) 

transpose  rvec,  store  as  tvec 

do  19  i=l,16 
ii  =  (i/4.1)+l 
alpha (ii)  =  O.dO 
do  19  j=l,2 

tvec(i,j)  =  rvec{j,i) 

19  continue 

normalize  eigenvector  matrix 

do  21  i^lflO 
ii  »  (i/4.1)+l 
do  21  j=l,2 

alpha(ii)  =  alpha(ii)  +  tvec(i,  j) »*2.d0 
21  continue 

do  23  i<.l,16 
ii  =  (i/4.1)+l 
do  23  j«l,2 

tnvec(i,j)  »  tvec {i,j) /dsqrt (alpha (ii)) 

23  continue 

write  (3,  •) 

write (3,*)  'normalized  eigenvectors  of  phi,  by  column' 
do  24  i=l,16 

write  (3, 7)  tnvec(i,l)  ,tnvec(i,2) 

24  continue 

7  format(lx,2(f20.13,lx)) 

compute  Poincare  exponents 
write  (3,*) 

write (3, »)  'Poincare  exponents' 
do  2100  i  =  1,4 
ww  •»  w(i) 

complex  log  of  eigenvalue  over  period 
xreal(i)  »  dlog(dsort(xww(l)*xww(l)  +  xww(2) *xww (2) ) ) 

•  /  period 

ximag(l)  "  dacan2(  xww(2),  xww(l)  )  /  period 
write (3, 5)  xreal (i) ,ximag(i) 

!100  continue 

5  format(3x,2(e20.13,lx)) 

create  input  file  for  flo.f  and  fho.f 

open (2, file“filnam2, status"' unknown' ) 

write (2,*)  xmu,'  ',xmua 
write (2,*)  period,'  ',npts 
write(2,*)  x(l,nxt),'  ',x(4,nxt) 
write  (2,*) 

write  normalized  eigenvector  parts  by  column 

—  if  complex  first  leave  alone 

—  if  real  first,  switch  order 

jt  -  0 

if  ((tnvcc(l,2) -cq.O.dO) -and. (tnvcc(2,2) -eq.O.dO))  jt=8 
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do  25  j"l,2 

do  25  i-l+jt,4+jt 

write  (2,*)  tnvec(i,j) 
25  continue 


c 

c 

c 

c 

c 


26 


do  26  i-9-jt,16-jt 
write(2,»)  tnvec(i,l) 
continue 

write  poincare  exponents  for  i  matrix 

—  if  complex  first  leave  alone 

—  if  real  first,  then  switch 


if  (ximagd) 

.eq.O.dO)  then 

write (2,*) 

write  (2,*) 

xreal  (3) 

write (2, ») 

ximag(4) 

write (2,*) 

O.dO 

write (2,*) 

O.dO 

write (2, •) 

ximag(3) 

write (2,*) 

xreal  (4) 

write  (2,*) 

O.dO 

write (2,*) 

O.dO 

write (2,*) 

O.dO 

write  (2,*) 

O.dO 

write (2, ») 

xreal  (1) 

write (2,*) 

ximag(2) 

write  (2,*) 

O.dO 

write  (2,*) 

O.dO 

write (2, •) 

ximagd) 

write  (2,*) 

xreal (2) 

else 

write  (2,  •) 
write  (2,  •) 

xreal d) 

write(2,  •) 

ximagd) 

write  (2, ») 

O.dO 

write  (2,*) 

O.dO 

write{2,*) 

ximagd) 

write  (2, ») 

xreal (2) 

write(2,*) 

O.dO 

write(2,*) 

O.dO 

write (2, •) 

O.dO 

write  (2,*) 

O.dO 

write (2, •) 

xreal  (3) 

write  (2,  •) 

xinag(4) 

write (2,*) 

O.dO 

write (2,*) 

O.dO 

write  (2,*) 

ximag(3) 

write  (2,  •) 

xreal (4) 

endif 

close (2) 
close (3) 


stop 

end 
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ccoccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  PROGRAM  FLOQUET/FOURIER 

<: 

c  PURPOSE:  Propagates  the  state  and  eigenvectors  of  a  periodic 

c  orbit  as  determined  by  program  period.  Every  value 

c  during  the  integration  is  saved,  and  then  converted 

c  into  a  fourler  series, 

c 

c  SUBROUTINES;  HAMING.F 

c  RHS2.F 

c  H.F 

c  FOURIER. F 

c 

Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
program  fl 

common  /lam/  xlambda(4) 
common  /data/  xmu,xmua,xj (4, 4) 

common  /ham/  t,x(20,4),£(20,4),err(20),nn,hh,mode 
c 

c  local  variables 

c 

implicit  double  precision  (a-h) 
implicit  inteoer  (1-n) 
implicit  double  precision  (o-z) 

dimension  xlambda  (4)  ,x(20,4)  ,£(20,4)  ,err  (20) ,  temp (100) ,  clc(51) 
dimension  xic  (2)  ,x  j  (4, 4)  ,x0  (16) ,  s  (4,  IfiO) ,  v  (4,  4, 100) ,  sk  (51) 
c 

c  read  input  data 

0 

read  (*,»)  xmu,xmua 
read  (*,*)  period, npts 
read  (»,»)  xic(l)  ,xlc(2) 

do  10  i=l,4 
do  10  j=l,4 

ii  -  (j-l)*4  +  i 
read  (*,»)  xO(ii) 

10  continue 

do  20  i-1,4 
do  20  1=1,4 

read  (*,»)  xj(j,i) 

20  continue 

hh  =  period/ (dble  (npts) ) 

write  (*,*)  'xmu  =  ',xmu,'  xmua  ■»  ',xmua 
write  (♦,*)  'orbit  period,  npts  ', period, npts 
write  (*,*)  'timestep  ',hh 

write  (*,*)  'initial  conditions  (ql,  pl=0,  q2=0,  p2) ' 
write  (*,*)  'ql  =  ',xic(l),'  p2  =  ',xic(2) 
write  (*,*) 

write  (*,*)  'f(O)' 
do  30  1=1,13,4 

write  (*,1)  x0(i),x0(i+l),x0(i+2)  ,x0(i+3) 

30  continue 

1  format (lx, £18.10, £18.10,  £18.10,  £18.10) 

write  (*,*) 
write  (*,*)  'xj(l,j)' 
do  40  i=l,4 

write  (*,1)  xj(i,l),xj(i,2),xj(i,3),xj(i,4) 

40  continue 

c 

c  set  up  initial  state 

c 

x(l,l)  =  xic(l) 
x(2,l)  =  O.dO 
x(3,l)  =  O.dO 
x(4,l)  »  xic(2) 
c 

c  initialize  £(0)  matrix 

do  160  i  -  1,16 

x(i+4,l)  =  x0(i) 

160  continue 

mode  =  1 
nn  =  20 
nxt  =  0 
t  =  O.dO 
c 

c  initialize  haming 
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499 


501 


502 

500 


call  hamlng(nxt) 

l£(nxt  .ne.  0)  go  to  499 

write  (*,*)  'failure  to  initialize' 
write  {»,*)  f{l,l),f(2,l) 
write  (»,*)  f(3,l),f(4,l) 
stop  99 
continue 

integration  loop 

l.do  500  i  =  1,100 
do  501  j  =  1,4 
s(j,i)  -  x(j,nxt) 
do  501  k  =  1,4 

v(j,k,i)  =  x(4*j+k,nxt) 
continue 
do  502  m  =  1,10 
call  hamlng(nxt) 
continue 
continue 

open  output  file 

open (2, file"' coef • fou' , status-' unknown' ) 


copy  eig  values/vectors  and  feed  to  fourier 

do  515  j=l,4 
do  510  i=l,100 
temp(i)  =  s(j,i) 

510  continue 

call  fourier (temp, ck,sk, 50) 
do  520  k=l,50 

write (2,*)  ck(k),sk(k) 

520  continue 

515  continue 


530 

535 

525 


do  525  1-1,4 
do  525  i=l,4 
do  530  k=l,100 

temp(k)  -  v(j,i,k) 
continue 

call  fourier (temp, ck,sk, 50) 
do  535  m-1,50 

write (2,*)  ck(m),sk(m) 
continue 


continue 


final  state  conditions 


write  (*,*) 

write (*,*)  'state  at  tf' 

writei*,*)  '  ql=' ,x  (l,nxt) , '  pl=',x(2,nxt) 
write(*,*)  '  q2=' ,x  (3,nxt) , '  p2=' ,x  (4,nxt) 

write  (*,*) 
write  (*,*)  'f(t)' 
do  600  i=5,17,4 

write  (*,1)  X (i, nxt) ,x (i+1, nxt) ,x (i+2,nxt) ,x (i+3,nxt) 
600  continue 


close  (2) 

stop 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  PROGRAM  FLOQUET/HAMILTONIAN 

c 

c  PURPOSE:  Calculate  period  coefficients  needed  In  the  new 

c  expanded  hamlltonlan,  from  the  third  order 

c  hamlltonlan  of  the  periodic  trajectory.  Compute 

c  after  each  Integration  step,  and  convert  the  result 

c  Into  a  fourler  series, 

c 

c  SUBROUTINES:  HAMING.F 

c  RHS2.F 

c  H.F 

c  FOURIER. F 

c 

Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


program  fh 

common  /lam/  x lambda (4) 
common  /data/  xmu,xmua,xj(4,4) 

common  /ham/  t,x(20,4) ,f (20,4), err(20),nn,hh, mode 
local  variables 


implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

dimension  xlambda  (4)  ,x  (20,4) ,  f  (20, 4) , err  (20)  ,xh3  (4,4,4) 
dimension  xic(2) ,  xj  (4, 4)  ,x0  (16)  ,xx  (4) ,  vl  (4) ,  v2  (4) 
dimension  tc(8),c(4,100)  ,c)c(100) ,  sic  (lOO) ,  temp  (100) 

read  input  data 

read  (»,*)  xmu,xmua 
read  (*,*)  period, npts 
read  (»,*)  xic(l)  ,xic(2) 


10 


transpose  col  to  row  to  fit  x(20) 

do  10  i=l,4 
do  10  j=l,4 

ii  -  (j-l)*4  +  i 
read  (*,*)  xO(ll) 
continue 


20 


do  20  i=l,4 
do  20  j=l,4 

read  (♦,»)  xj(j,i) 
continue 


hh  =  period/ (dble (npts) ) 
output  inputs 

write  (*,*)  'xmu  »  ',xmu,'  xmua  »  ',xmua 
write  (*,*)  'orbit  period,  npts  ', period, npts 
write  (*,*)  'timestep  ',hh 

write  (*,*)  'initial  conditions  (ql,  pl“0,  q2=0,  p2)' 
write  (*,*)  'ql  =  ',xic(l),'  p2  =  ',xic(2) 
write  (*,*) 

write  (»,*)  'f(O)' 
do  30  1=1,13,4 

write  (*,1)  x0(i),x0(i+l),x0(i+2),x0(i+3) 

30  continue 

1  format (lx, 4 (fl8. 10)) 

write  (*,*) 

write  (*,*)  'xj(i,j)' 

do  40  i=l,4 

write  (*,1)  xj(i,l),xj(i,2),xj(i,3),xj(i,4) 

40  continue 

set  up  initial  state 

x(l,l)  =  xic(l) 

x(2,l)  =  O.dO 

x(3,l)  =  O.dO 

x(4,l)  =  xlc(2) 

initialize  f(0)  matrix 

do  160  i  =  1,16 
x(i+4,l)  »  x0(i) 

160  continue 
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mode  ■=  1 
nn  “  20 
nxt  “  0 
t  “  O.dO 

Initialize  hamlng 

call  hamlng (nxt) 

if (nxt  .ne.  0)  go  to  499 

write  (*,*)  'failure  to  initialize' 
write  (*,*)  f(l,l),f(2,l) 
write  (*,*)  f(3,l),f(4,l) 
stop  99 

499  continue 

begin  integration  loop 

do  500  i  =  1,100 
do  510  j  =  1,4 
xx(j)  -  x(j,nxt) 
vl(j)  »  x(4*j+l,nxt) 
v2(3)  "  x(4*j+2,nxt) 

510  continue 

compute  third  order  hamiltonlan  tensor 

do  520  j  -  1,4 
do  520  k  =  1,4 
do  520  m  =  1,4 

xh3(j,k,m)  -  h(xx,3,  j,k,m,0,0) 

520  continue 

do  525  j=l,8 
tc(j)  =  O.dO 
525  continue 

compute  periodic  coefficients 

do  530  j=l,4 
do  530  k=l,4 
do  530  m=l,4 

tc(l)  =  tc(l)  +  xh3(j,k,m)  »  vl(j)  *  vl  (k)  *  vl  (m) 

tc(2)  =  tc(2)  +  xh3(3»k,m)  *  vl(3)  *  v2(k)  *  vl  (m) 

tc(3)  »•  tc(J)  +  xh3(3»k,m)  *  v2(3)  *  vl  (k)  *  vl  (m) 

to(4)  =  tc(4)  +  xh3(3»k,m)  *  vl(3)  »  vl(k)  *  v2  (m) 

tc(5)  =  tc(5)  +  xh3(j,k,m)  *  v2(j)  •  v2(k)  *  vl  (m) 

to(6)  -  tc(6)  +  xh3(j,k,m)  »  vl(j)  *  v2(k)  »  v2(m) 

tc(7)  -  tc(7)  +  xh3(j,k,m)  »  v2(3)  *  vl  (k)  *  v2  (m) 

tc(8)  =  tc(8)  +  xh3(j,k,m)  *  v2(j)  *  v2(k)  *  v2(m) 
530  continue 

c(l,i)  =  tc(l)/6.d0 

c(2,i)  =  (tc(2)+tc(3)+tc(4))/6.d0 

c(3,i)  =  (tc(5)+tc(6)+tc(7))/6.d0 

c(4,i)  =  tc(8)/6.d0 

do  550  m  °  1,10 

call  hamlng (nxt) 

550  continue 

500  continue 

compute  fourier  coefficients  from  periodic  ones 

open (2, file='coef .ham' , status"' unknown' ) 

do  570  1=1,4 
do  580  j=l,100 
temp(3)  "  c(i,j) 

580  continue 

call  fourier (temp, ck, sk, 50) 

do  590  k=l,50 

write(2,*)  ck(k),sk(k) 

590  continue 

write (2, *) 

570  continue 

final  state  conditions 
write  (»,*) 

write (*,*)  'state  at  tf' 

write(»,*)  '  ql"',x(l,nxt),'  pl»' ,x(2,nxt) 
writel*,*)  '  q2"',x(3,nxt),'  p2»' , x (4, nxt) 
write  (*,♦) 
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write  (*,*)  'f(t)' 
do  600  i=S,17,4 

write  (*,1)  x(i,nxt)  ,x(i+l,nxt)  ,x(i+2,nxt)  ,x(i+3,nxt) 
continue 

close  (2) 

stop 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  PROGRAM  EXACT 

c 

c  PURPOSE:  Integrates  a  nearly  periodic  trajectory,  subtracts 

c  the  periodic  reference,  and  transforms  the  result 

c  into  modal  variables,  and  creates  a  plotfile 

c  bl  vs  b2. 

c 

c  SUBROUTINES:  HAMING.F 

c  RHSl.F 

c  H.F 

c 

Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

program  check 

problem  commons 

common  /data/  xmu,xmua 
common  /lam/  xlambda(4) 

common  /ham/  t,x(20,4),f (20,4),err(20),nn,hh,mode 

variable  declarations 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 

Implicit  double  precision  {o-z) 
character*10  fi Inaml, fllnam2, filnaml 

dimension  xlambda (4) ,x (20, 4) , f (20, 4) , err (20) 
dimension  slnn (50) ,  coss (SO) ,  ck (20, 50) , sk (20, 50) , cf (20) 
dimension  cc(4,  4)  ,xx(4)  ,dx(4)  ,xxx(50) 

input  data 

read(*,*)  xmu,xmua 
tead(*,*)  period, npts 
hh  =  period/ (dble (npts) ) 
read(*,*)  xnot,ynot 
read(*,*)  xjac, syn 
read(*,*)  w,trip 
read(*,*)  filnaml 
read(*,*)  filnam2 
read(*,»)  filnam3 


do  100  i=l,20 
do  100  j-l,S0 

read(*,»)  ck(i, j),sk(i, j) 

100  continue 

mode  =  0 
nn  -  4 
nxt  =  0 
t  =  O.dO 

pi  =  dacos(-l.dO) 
wO  “  2.d0*pi/period 

get  ql,pl,q2,p2  for  given  xO,yO,  and  jacobian 

ql  =  xnot  +  xrou 
q2  =  ynot 

xham  =  (xrou*xmua-x jac) /2.d0 
rl  =  ((ql-xmu)*»2.d0  +  q2**2.d0) »» .5d0 
r2  =  ( (ql+xmua) **2.d0  +  q2**2.d0) ‘‘.SdO 
d  =  xham  +  xmua/rl  +  xmu/r2 
g  =  q2/(ql-xmu) 
a  “  g*g  +  i.dO 

b  =  -2.d0* (g*g*xmu  +  g*q2  +  ql) 
c  =  (g*xmu) **2.d0  +  2.d0*g*q2*xmu  -  2.d0*d 
p2  =  (-b+syn*(b*L-4.dO»a*c)»».5dO)/(2.dO*a) 
pi  =  g*(xmu-p2) 

write (*,*)  'initial  conditions  (Jefferys) ' 
write(*,*)  'xO»',xnot,'  yO“',ynot 
write ♦) 

write (*!*)  'initial  conditions  (Szebeheiy) ' 
write  (*,*)  'ql=',ql,'  pl>»',pl 
write  (*,*)  'q2=',q2,'  p2-',p2 
write  (*,  *) 

write(*,*)  '  jacobian"' , xjac 
write  {*,*)  'period"' , period 
write  (*,  •) 

write  (*,*)  'initial  conditions  (Modal)' 
c 
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initial  conditions 


x(l,l) 

x(2,l) 

x(3,l) 

x(4,l) 


p2 


initialize  haming 
call  haming (nxt) 
turn  off  second  EOM  eval 


nxt  =  -nxt 

if  (nxt  .ne.  0)  go  to  499 
stop  99 

499  continue 

open  output  files 

open (2, file»filnaml, status”' unknown' ) 
open (3, file=filnam2, status”' unknown' ) 
open (4,  file=filnam3, status”' unknown' ) 

integration  loop 

do  500  i  =  0,npts»trlp 

if  (mod(l,20)  .eq.O)  then 


compute  sin  (n*theta) ,  cos {n*theta) ,  n=l  to  50 

cossd)  “  dcos(wO*t) 

sinn{l)  =  dsin(wO*t) 

coss(2)  =  2.d0»coss(l)*coss(l)  -  l.dO 

slnn(2)  =  2.d0*sinn(l)*coss(l) 

do  200  j=3,S0 

coss(j)  ”  2-d0*coss(j-l)*coss(l)  -  coss(j-2) 
sinn(j)  =  2.d0»sinn(j-l)»coss(l)  -  sinn(3-2) 
200  continue 

reassemble  periodic  traj  and  eigenvector  matrix 

do  300  k»l,20 
cf(k)  ”  ck(k,l) 

do  300  j=l,49 

cf(k)  ”  cf(k)  +  ck(k,  j+l)*coss(j) 

•  +  sk(k, j+l)»sinn(j) 

300  continue 

write  (*, *) 
do  301  j  =  1,20 
write{*,»)  cf(j) 

301  continue 


compute  generalized  eigenvector  (grad  of  hamiltonian) 

do  320  j”l,4 
xx(j)  =  cf(j) 
continue 

temp  ”  O.dO 
do  330  j=l,4 

dx(j)  =  x(j,  iabs(nxt) )  -  cf(j) 
cf(I6+j)  ”  h(xx,l,i,0,0,0,0) 
temp  =  temp  +  cf  (16+ j) ‘cf  (16+ j) 
continue 

write(*,*) 
do  331  j=l,4 
write(+,*)  cf{16+j) 
continue 


do  335  j”l,4 

cf(16+j)  =  cf  (16+j)/(dsqrt(torop) ) 
continue 

writc(*,») 
do  336  j-1,20 
write(*,*)  cf(j) 
continue 


320 


330 

Un 


335 


c 

c 

c 

c  336 
c 
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place  eigenvectors  in  4x4  matrix  (for  inversion) 

do  340  j-1,4 
do  340  k=1.4 
jj  =  4*j+)c 
co(k,j)  =  cf(jj) 

340  continue 

c 

c  invert  delta  x  =  eigenvector  matrix  *  b 

c 

idig  -  0 

call  leqt2£ (cc, 1, 4, 4,dx, idlg,xxx, ier) 

if  (i.eq.O.dO)  then 

write{*,*)  'bl0-',dx(l),'  b20-'.dx(2) 
endif 

write(2,»)  dx(l).dx(2) 
write(3,*)  t,dx(3) 
write(4,*)  t,dx(4) 

endif 

call  haming(nxt) 

500  continue 

close (2) 

stop 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  PROGRAM  EXPANDED 
c 

c  PURPOSE:  Using  the  periodic  coefficients  made  by  the  program 
c  flo^et/hamlltonlan,  the  eom  for  the  truncated 

c  hamlltonlan  case  are  integrated.  A  plotflle 

c  matching  bl  vs  b2  is  created, 

c 

c  SUBROUTINES:  NAMING. F 

c  RHS3.F 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
program  )c3 

common  /data/  wO, w, ck (4,50) , sk (4, 50) 
common  /ham/  t,x(20,4),f {20,4),err(20),nn,hh,mode 
c 

c  local  variables 

c 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
character*10  filnam 

dimension  x  (20,4) ,  f  (20, 4)  ,etr  (20) 
dimension  ck  (4, 50) ,  sk  (4,50) 

read  input  data 

read  (*,*)  period, npts 
read  (*,*)  bl0,b20 
read  (»,*)  w,trip 
read(»,*)  filnam 


read  fourier  coefficients 

do  20  i=l,4 
do  20  j=l,50 

read  (*,»)  ck(i,  j)  ,sk(i,  j) 

20  continue 

hh  =  period/ (dble  (npts) ) 

output  inputs 

write  (•,*)  'orbit  period,  npts  ', period, npts 
write  (*,*)  'timestep  ',hh 
write  (*,»)  'initial  conditions  (modal)' 
write  (*,»)  'bl0=',bl0,'  b20=',b20 

set  up  initial  state 

x(l,l)  »  blO 
x(2,l)  =  b20 

mode  “  0 
nn  =  4 
nxt  =  0 
t  «  O.dO 

pi  =  dacos(-l.dO) 
wO  =  2.d0*pi/period 

initialize  haming 

call  haming (nxt) 

if  (nxt  .ne.  0)  go  to  499 

write  (*,*)  'failure  to  initialize' 
write  (»,*)  f (1,1), f (2,1) 
write  (»,*)  f(3,l),f(4,l) 
stop  99 

499  continue 

begin  integration  loop 

open (2, filesfilnam, status-' unknown' ) 

do  500  i  =  l,npts*trip 
call  haming (nxt) 
if  (mod(i,l00) .nc.O)  go  to  500 
writc(2,»)  X  (l,nxt)  ,x(2,nxt) 

500  continue 
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close  (2) 

stop 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  SUBROUTINE  NAMING 

c 

c  PURPOSE:  Hatnlng  is  an  ordinary  differential  equations 

c  inteqrator.  It  is  a  fourth  order 

c  predictor-corrector  algorithm  that  carries  the  last 

c  four  values  of  the  state  vector,  extrapolates 

c  them  to  obtain  the  next  value  (the  prediction  part), 

c  and  then  corrects  the  extrapolated  value  to  find  a 

c  new  value  for  the  state  vector.  Nxt  specifies  which 

c  of  the  4  values  of  the  state  vector  is  the  "next" 

c  one.  Nxt  is  updated  by  haming  automatically,  and  is 

c  zero  on  the  first  call.  The  user  supplies  an 

c  external  routine  rhsCnxt),  which  evaluates  the 

c  equations  of  motion 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  haming (nxt) 

common  /ham/  x,y(20, 4),f (2C,4),errest(20),n,h,mode 
implicit  double  precision  (a-z) 
dimension  y  (20, 4 ) , f (20, 4) ,errest (20) 
integer  l,l,nxt,n,npl,nml,nm2,npo,isw,  jsw,made 
c 

c  X  is  the  Independent  variable  (  time  ) 

c  y(6, 4)  is  the  state  vector-  4  copies  of  it,  with  nxt 

c  pointing  at  the  next  one 

c  f(6,4)  are  the  equations  of  motion,  again  four  copies 

c  a  call  to  rhs(nxt)  updates  an  entry  in  f 

c  errest  is  an  estimate  of  the  truncation  error  -  normally  not 
c  used 

c  n  is  the  number  of  equations  being  integrated  -  6  or  42  here 

c  h  is  the  time  step 

c  mode  is  0  for  just  EOM,  1  for  both  EOM  and  EOV 

c 

tol  ■>  0.0000000001 

if  (nxt)  190,10,200 

10  xo  ■»  X 

hh  «  h/2.0d+00 
call  rhs(l) 
do  40  1  =  2,4 
X  =  X  +  hh 
do  20  i  =  l,n 

20  y(i,l)  -  y(i,l-l)  +  hh*f(i,l-l) 
call  rhs(i) 

X  »  X  +  hh 
do  30  1  a  l,n 

30  y(i,l)  -  Y(i,l-1)  +  h»f{i,l) 

40  call  rhs(l) 
jsw  -  -10 
SO  isw  -  1 

do  120  1  >1,0 

hh  -  y(i,l)  +  h*(  9.0d+00»f (i.l)  *  19.0d+00* f (i, 2) 

1  -  5.0d+00»f (i,3)  +  f(i,4)  )  /  24.0d+00 

if(  dabs(  hh  -  y(i,2))  .It.  tol  )  go  to  70 
isw  =  0 

70  y(i,2)  =  hh 

hh  -  y(i,l)  +  h»(  f(i,l)  +  4.0dt00»f  (i,2)  +  f  (i,  3) ) /3.0d+00 
if(  dabs(  hh-y(i,3))  .It.  tol  )  go  to  90 
isw  B  0 

90  y(i,3)  =  hh 

hh  -  y(i,l)  +  h»(  3.0d+00»f (i,l)  ♦  9.0d+00*f (i.2)  + 

1  9.0df00»f(i,3)  +  3.0d+00»f(i,4)  )  /  8.0d+00 

if(  dabs  (hh-y  U,  4) )  .It.  tol  )  go  to  110 
isw  »  0 

110  y(i,4)  =  hh 
120  continue 
X  "  xo 

do  130  1-2,4 
X  -  X  +  h 
130  call  ths(l) 

if(isw)  140,140,150 
140  jsw  -  jsw  ♦  1 

if(jsw)  50,280,280 
ISO  X  -  xo 
isw  -  1 
jsw  »  1 

do  160  i  -  i,n 
160  crrcst(i)  -  0.0 
nxt  -  1 
go  to  280 
190  jsw  =  2 
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nxt  -  labs(nxt) 
c 

c  this  is  hamings  normal  propagation  loop  - 
c 

200  X  =  X  +  h 

npl  “  mod (nxt, 4)  +  1 
go  to  (210,230) , isw 
c  permute  the  index  nxt  modulo  4 
210  go  to  (270, 270, 270, 220), nxt 
220  isw  =  2 

230  nm2  mod  (npl,  4)  +  1 
nml  “  mod(nm2,4)  +  1 
npo  “  mod (nml, 4)  +  1 

this  is  the  predictor  part 

do  240  i  -  l,n 

f(i,nm2)  =  y(i,npl)  +  4.0d+00»h*(  2.0d+00»f (i,npo)  -  f(i,nml) 
1  +  2.0d+00»f (i,nm2)  )  /  3.0d+00 

240  >(i,npl)  =  f(i,nm2)  -  0.92S619835*errest  (i) 


c  now  the  corrector  -  fix  up  the  extrapolated  state 
c  based  on  the  better  value  of  the  equations  of  motion 
c 

call  rhs(npl) 
do  2S0  i  »  l,n 

y(i,npl)  =  (  9.0d+00*y (i,npo)  -  y(i,nm2)  +  3.0d+00*h*(  f(i,npl) 
1  +  2.0dt00*f (i,npo)  -  f(i,nml)  )  )  /  8.0d+00 

errest(i)  =  f(i,nm2)  -  y(i,npl) 

250  y(i,npl)  =  y(i,npl)  +  0.0743801653  *  errest(i) 
go  to  (260,270) , jsw 
260  call  r)is(npl) 

270  nxt  =  npl 
280  return 
end 
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000  000  OOO  OOO  O  000  OOO  O  O  000 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  SUBROUTINE  RHSl 
c 

c  PURPOSE:  Computes  the  right-hand  side  of  the  equations  of 
c  motion  and  variation,  in  the  restricted  three-body 

c  problem, 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  rhs (k) 

canonical  EOM  and  EOV,  4th  order  system 

common  /lam/  xlambda(4) 
common  /data/  xmu,  xmua 

common  /ham/  t,  x{20,4) ,f (20,4) ,err{20) ,nn,hh,mode 

double  precision  t,x, f,err,hh 
double  precisian  xlambda, xmu, xmua 

double  precision  h,xx (4) ,z (4,4) ,grdh (4, 4) , temp(4, 4) 


data 

z/  O.dO, 

-l.dO. 

O.dO, 

O.dO, 

1 

l.dO, 

O.dO, 

O.dO, 

O.dO, 

2 

O.dO. 

O.dO, 

O.dO, 

-l.dO, 

3 

O.dO, 

O.dO, 

l.dO, 

O.dO/ 

extract  state 

do  10  i  =  1,4 
xx{i)  =  x(i,k) 

10  continue 

equations  of  motion 

f{l,k)  «  h(xx,l, 2,0, 0,0,0) 
f(2,k)  -  -h(xx,i.l,0,0,0,0) 
f(3,k)  •>  h(xx,l,4,0,0,0,0) 
f(4,k)  =  -h  (XX,  1,3, 0,0, 0,0) 

if (mode  .eq.  0)  return 

calculate  order  2  gradient  matrix 

do  20  1  •>  1,4 

do  20  j  -  1,4 

grdh(i,j)  =  h(xx,2,i,  j,0,0,0) 

20  continue 

matrix  mpy  by  z 

do  30  i  »  1,4 

do  30  ii  °  1,4 
ij  =  4*i+ii 
temp(i,ii)  =  O.dO 
do  3()  j  =  1,4 

temp(i,ii)  «•  temp(i,ii)  +  z (i,  j)  »grdh ( j,  ii) 

30  continue 

calculate  A  phi 

do  35  i  =  1,4 

do  35  ii  ■  1,4 
ij  «  4»i+ii 
f(ii,k)  =  O.dO 
do  35  j  =  1,4 

f(ij,k)  -  f(ij,k)  +  tcmp(i,  j)»x(4*  j+ii,k) 

35  continue 

check  if  propagating  eigenvectors  or  phi 

if  (mode  .cq.  1)  return 
c 

c  get  eigenvector  com 

do  40  j  -  2,4,2 

do  40  ii  -  1,4 

ij  «>  4»(j-i)  +  II 

f(ij,k)  »  f(ij,k)  +  xlambda (j)»x(ij+4,k) 
f(i3+4,k)  -  f{ij+4,k)  -  xla-bda(j)*x(i j,k) 

40  continue 
return 
end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  SUBROUTINE  RHS2 

c 

c  PURPOSE:  Creates  right-hand  side  of  the  eos  for  the 
c  restricted  three-body  problem,  and  for  the 

c  eigenvector  equation, 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine  rhs(k) 
c 

c  canonical  EOM  ar.d  EOV,  4th  order  system 

c 

common  /lam/  xlainbda(4) 
common  /data/  xnu,xnua,xj (4,4) 

common  /ham/  t,x(20,4) , f  (20, 4) ,err (20) ,nn,hh,made 

double  precision  x  (20, 4) ,  f  (20, 4)  ,err  (20) ,  fl  (20, 4) ,  f2  (20,  4) 
double  precision  xlambda,xmu,xmua,t,hh 

double  precision  h,xx(4),z(4,4),grdh(4,4),temp(4,4),xj(4,4) 


data 

2/  O.dC, 

-l.dO, 

O.dO, 

O.dO, 

1 

l.dO, 

O.dO, 

O.dO, 

O.dO, 

2 

O.dO, 

O.dO, 

O.dO, 

-l.dO, 

3 

O.dO, 

O.dO, 

l.dO, 

O.dO/ 

extract  state 

do  10  i  =  1,4 
xx(i)  »  x(i,)c) 

10  continue 

equations  of  motion 

f(l,k)  »  h(xx,l,2,0,0,0,0) 
f(2,k)  =  -h(xx,l, 1,0, 0,0,0) 
f(3,k)  -  h(xx,l, 4,0, 0,0,0) 
f(4,k)  =  -h(xx,l,3,0,0,0,0) 

if (mode  .eq.  0)  return 

calculate  order  2  gradient  matrix 

do  20  i  -  1,4 

do  20  j  -  1,4 

grdh(i,j)  =  h(xx,2,i,  j, 0,0,0) 

20  continue 

matrix  mpy  by  z 

do  30  i  •>  1,4 

do  30  ii  o  1,4 
ij  -  4*i+ii 
cemp(i,ll)  =  O.dO 
do  jO  j  -  1,4 

tempu.ii)  “  tecp(i,ii)  *  z(i,  j)  •grdh(j,ii) 

30  continue 

calculate  A  phi  and  phi  J 

do  35  i  -  1,4 

do  35  ii  -  1,4 
i1  -  4*i+ii 
fl(ij,k)  -  O.dO 
f2(i1,k)  -  O.dO 
do  35  i  -  1,4 

f2(ij,k)  -  f2(ij,k)  ♦  x(4*i+j,k)»xj(j,ii) 
fl(ij,k)  »  fl(ij,k)  +  tcmp(i,  j)»x(4»j+ii,k) 


35  continue 


do  36  i  -  5,20 

f(i,k)  -  a(i,k)  -  f2(i,k) 
36  continue 


return 

end 
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OOO  OOft  ooo  ooo  oooooo 


ccccccccccccccccccccccccccecccccccccccccccccccccccccccccccoccccc 

SUBROUTINE  RHS3 

PURPOSE;  Calculate  rhs  for  nearly-periodic  eom,  using 
expanded  hamlltonlan. 


subroutine  rhs (k) 

canonical  EOM  and  EOV,  4th  order  system 

common  /data/  w0,w,ck(4,50),sk(4,50) 

common  /ham/  t,x (20, 4) , f (20, 4) ,err (20) ,nn,hh,mode 

double  precision  t,x (20, 4) , f (20,4) ,err (20) ,hh, sinn (50) 
double  precision  ck (4, 50) , sk (4, 50) , c (4) , coss (50) , wO, w 
doul^le  precision  bl,b2 


generate  sin(l  to  SO  *  wO)  and  cos(l  to  50  *  i;0) 
coss(l)  “  dcos(w0*t) 

coss(2)  "  2-d0*co33(l)*coss(l)  -  l.dO 

slnn(l)  ■»  dsin(w0*t) 

slnn(2)  “  2.d0*slnn(l) *coss(l) 

do  100  i=3,50 

coss(i)  =  2.d0*coss (i-1) *coss (1)  -  coss(i-2) 
slnn(i)  =  2.d0*sinn (i-l) *coss (1)  -  sinn(i-2) 
100  continue 


reconstruct  periodic  function  from  coefficients 

do  200  i=l,4 
c(i)  =  ck(i,l) 
do  200  j=l,49 

c(i)  =  c(i)  +  ck(i,  j+l)’^coss(j)  +  sk(i,  j+l)*sinn(j) 
200  continue 

bl  =  x(l,k) 
b2  =  x(2,k) 

calculate  bl  dot  and  b2  dot 

f(l,k)  =  w*b2  -  bl*bl*c(2)  + 

*  2.d0*bl*b2*o(3)  -  3.d0*b2*b2*o (4) 

£(2,k)  =  -w*bl  -  b2*b2*c(3)  + 

*  2.d0*bl*b2*c(2)  -  3.d0*bl*bl*c(l) 

return 

end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c  FUNCTION  H 
c 

c  PURPOSE;  Computes  desired  order  tensor  of  the  restricted 
c  three-body  hamlltonlan. 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


function  h(x,lord,l, 
c 

c  state  vector  x  =  (  ql,  pi,  q2,  p2  ) 
c 

common  /data/  xmu,xmua 
Implicit  double  precision  (a-z) 
dimension  x(4) 

Integer  lord, jord,l, j,k 
c 

c  preliminaries 
c 

qa  ■>  x(l)  -  xmu 
qb  =  x(l)  +  xmua 
rl  =  dsqrt  (qa*qa  +  x(3)*x(3)) 
r2  =  dsqrt (qb*qb  +  x(3)*x(3)) 
c 

c  branch  on  order 
c 

jord  =  lord  +  1 

go  to  (1,  1000,  2000,  3000) , jord 
c 

c  **  Order  Zero  »♦ 

c 

1  continue 

h  -  0.5d0*(x(2)*x(2)  +  x(4)*x(4))  +  x(3)*x(2)  -  x(l)*x(4) 
1  -  xmua/rl  -  xmu/r2 

return 

1000  continue 
c 

c  **  Order  One  ** 

rl3  =  rl**3.d0 
r23  =  r2**3.d0 

go  to  (1001,  1002,  1003,  1004),  1 

1001  h  =  -x(4)  +  xmua*qa/rl3  +  xmu*qb/r23 
return 

1002  h  =  x(2)  +  x(3) 
return 

1003  h  =  x<2)  +  xmua*x(3)/rl3  +  xmu*x(3)/r23 
return 

1004  h  =  x(4)  -  x(l) 
return 


2000  continue 
c 

c  *»  Order  Two  ** 

rl3  =  rl**3.d0 
r23  =  r2»*3.d0 
rl5  =  rl**5.d0 
r25  =  r2**5.d0 

go  to  (2001,  2002,  2003,  2004),! 

2001  go  to  (2011,  2012,  2013,  2014),  j 

2002  go  to  (2021,  2022,  2023,  2024), j 

2003  go  to  (2031,  2032,  2033,  2034) , j 

2004  go  to  (2041,  2042,  2043,  2044), j 

2011  h  =  xmua/rl3  +  xmu/r23  -3.d0*xmua»qa»qa/rl5 

1  -3.d0*xmu*qb*qb/t25 

return 

2012  h  =  O.dO 
return 

2013  h  “  -3 .d0*xmua*qa*x (3) /rl5  -  3.d0*xmu*qb*x (3) /r25 
return 

2014  h  =  -l.dO 
return 

2021  h  =  O.dO 
return 

2022  h  =  l.dO 
return 

2023  h  =  l.dO 
return 

2024  h  =  O.dO 
return 

2031  go  to  2013 
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2032  h  =  l.dO 
return 

2033  h  =  -3,d0*xmua*x(3)*x(3)/rl5  -  3.d0*xmu»x{3)*x(3)/r25 

1  +  xmua/rl3  +  xmu/r23 

return 

2034  h  =  O.dO 
return 

2041  h  =  -l.dO 
return 

2042  h  =  O.dO 
return 

2043  h  =  O.dO 
return 

2044  h  =  l.dO 
return 


3000  continue 
c 

c  *•  Ort  er  Three 

rl5  =  rl**S.dO 
r25  =  r2**5.d0 
rl7  =  rl**7.d0 
r27  =  r2**7.d0 


go  to  (30001,  30002,  30003,  30004), i 

30001  go  to  (30110,  30120,  30130,  30140) , j 

30002  go  to  (30210,  30220,  30230,  30240),  j 

30003  go  to  (30310,  30320,  30330,  30340), j 

30004  go  to  (30410,  30420,  30430,  30440) , j 

c  note  matrix  Is  quite  sparse  now . 

30110  go  to  (30111,  30112,  30113,  30114), )c 

30130  go  to  (30131,  30132,  30133,  30134), Ic 

30310  go  to  (30311,  30312,  30313,  30314), )c 

30330  go  to  (30331,  30332,  30333,  30334), )c 


30111  h  =  -9.d0*xmua*qa/rlS  -  9.d0*xmu*qb/r25 

1  +  15.d0*xmua*qa*qa*qa/rl7  +  15.a0*xmu*qb*qb*qb/r27 

return 

30112  h  =  O.dO 
return 

30113  h  =  -3.d0*xmua*x(3)/rl5  -  3.d0*xmu*x(3) /r25 

1  +  lS.dO*xmua*qa*qa*x(3) /rl7  +  lS.d0*xmu*qb*qb*x (3) /r27 

return 

30114  h  =  O.dO 
return 

30120  h  =»  O.dO 
return 

30131  go  to  30113 

30132  h  =  O.dO 
return 

30133  h  =  -3.d0*xmua»qa/rl5  -  3.d0*xmu»qb/r25 

1  +  15.d0*xmua*qa*x(3)*x(3)/rl7  +  15.d0»xmu*qb*x (3) *x (3) /r27 

return 

30134  h  =  O.dO 
return 

30140  h  =  O.dO 
return 

30210  h  =  O.dO 
return 

30220  h  =  O.dO 
return 

30230  h  ■=  O.dO 
return 

30240  h  =  O.dO 
return 

30311  go  to  30113 

30312  h  =  O.dO 
return 

30313  go  to  30133 

30314  h  =  O.dO 
return 

30320  h  =  O.dO 
return 

30331  go  to  30133 

30332  h  =>  O.dO 
return 

30333  h  =  -9.d0*xmua*x(3)/rl5  -  9.d0*xmu*x (3) /r25 

1  +  IS.dO* (xmua/rl7  +  xmu/r27) *x (3) ‘x (3) *x (3) 

return 

30334  h  =  O.dO 
return 

30340  h  "  O.dO 
return 

30410  h  =  O.dO 
return 
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30420  h  =  O.dO 
return 

30430  h  »  O.dO 
return 

30440  h  =  O.dO 

return 

end 


QOQ  000  OOOOOO 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  SUBROUTINE  FOURIER 

c 

c  PURPOSE:  harmonic  analysis  of  2n  values  of  function  F 
c  evenly  spaced  at  Interval  2pl/2n,  starting  with 

c  zero.  Into  n+1  cosine  coefficients  ck  and  n-1  sine 

c  coefficients  sk. 

c 

c  ref  Brouwer  and  Clemence,  p  109 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine  fourler(F,ck,sk,  n) 

double  precision  F(2) ,ck(2) , sk(2) , twopl,alpha 

twopl  =  2.d0*3.141592653589d0 
alpha  =  twopl/dble (2*n) 
n2ml  “  2*n-l 

order  k  loop 

do  500  k  =  0,n 

cosine  sum 

ck(k+l)  -  O.dO 

do  200  j  =  0,n2ml 

ck(k+l)  =  ck{k+l)  +  F(j+1)  *  dcos(  dble (k» j) *alpha) 
200  continue 

ok(k+l)  =  ck (k+i) /dble  (n) 

sine  sum 

lf(k  .eq.  0)  go  to  500 
lf(k  .eq.  n)  go  to  500 
sk(k+l)  =  O.dO 
do  400  j  l,n2ml 

sk(k+l)  -  sk(kH)  +  F(j+1)  *  dsln(  dble  (k*j) ‘alpha) 
400  continue 

sk(k+l)  =  sk(k+l)/dble(n) 

500  continue 

correct  first  and  last  cosine  coefficient 

ck(l)  =  0.5d0»ck{l) 
ck(n+l)  =  0.5d0*ck(n+l) 

return 
end 
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